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ABSTRACT 


To  successfully  combine  information  from  distributed  radar  sensors,  it  is  essential  that  each 
sensor  be  correctly  referenced  in  a  global  coordinate  system.  If  there  are  biases  in  the  reported  position  of 
these  sensors,  the  reported  target  position  will  also  be  biased  and  the  ensuing  global  estimate  of  the  target 
position  will  be  degraded.  Furthermore,  any  biases  in  range  or  azimuth  measurements  of  these  sensors 
will  likewise  be  reflected  in  the  degradation  of  global  estimate  of  target  position.  Registration  is  the 
process  of  ensuring  that  these  errors  do  not  result  in  the  creation  of  an  additional  redundant  target  when 
only  a  single  target  exists. 

The  objective  of  this  thesis  is  to  create  a  model  for  analyzing  the  impact  of  these  biases 
quantitatively.  The  model  consists  of  modules  which  perform  the  required  coordinate  conversion, 
tracking,  and  data  correlation.  The  target  tracks  are  provided  by  a  standard  Kalman  filter  assuming  a 
constant  velocity  model.  The  measurements,  state  estimates,  and  covariance  matrices  obtained  from  the 
Kalman  filter  are  combined  to  form  a  Chi-squared  correlation  gate.  With  this  model,  the  bounds  on 
position,  range,  and  azimuth  biases  are  determined  individually  and  cumulatively.  The  simulated  results 
compare  favorably  with  the  theoretically  determined  bounds.  An  additional  benefit  of  this  model  is  that 
the  spatial  dependence  of  the  biases  may  be  obtained. 
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THE  EFFECT  OF  REGISTRATION  ERRORS 
ON  TRACKING  IN  A 
NETWORKED  RADAR  SYSTEM 


1.  Introduction 


Air  defense  and  air  traffic  control  systems  depend  on  a  surveillance  subsystem  to  provide  an 
overall  picture  of  the  air  situation  in  order  to  make  decisions.  To  maintain  an  accurate,  complete,  and 
current  air  picture,  the  surveillance  subsystem,  in  turn,  depends  on  combinations  of  networked  sensors  to 
provide  the  raw  data  from  which  the  picture  is  constructed. 

The  general  registration  problem  arises  whenever  it  is  desired  to  combine  information  from  two 
or  more  sensors  into  a  single  “system  level”  surveillance  picture.  The  most  important  attribute  of  a 
“good”  surveillance  picture  is  that  it  contains  exactly  one  track  for  each  object  detected  by  at  least  one 
sensor  in  the  system.  The  fundamental  problem  in  sensor  networking,  therefore,  is  to  determine  whether 
the  data  reported  by  two  or  more  sensors  represent  a  common  object  or  two  (or  more)  distinct  objects. 

There  are  often  systematic  errors  which  effectively  introduce  biases  into  the  track  estimation 
process.  Therefore,  failure  to  register  adequately  in  a  multiple  radar  system  can  result  in  varying  degrees 
of  performance  degradation,  depending  on  the  magnitude  of  the  biases  with  respect  to  the  random 
measurement  errors  and  the  track  correlation  gates.  The  level  of  degradation  ranges,  at  worst,  from  the 
formation  of  multiple  redundant  tracks  for  a  single  aircraft  to  reduced  track  accuracy  and  stability  when 
the  bias  is  relatively  small.  In  between,  the  benefits  of  a  multiple  radar  system  can  be  negated  and  the 
system,  in  effect,  can  be  reduced  to  a  single  radar  tracking  system[Bar90], 
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It  should  be  noted  that  this  topic  is  really  a  special  case  of  a  more  general  sensor  fusion  problem 
which  may  employ  many  diverse  sensors.  Data  from  different  sensors  and  types  of  sensors  are  combined 
using  techniques  drawn  from  several  disciplines,  such  as  signal  processing,  statistics,  artificial 
intelligence,  pattern  recognition,  cognitive  psychology,  and  information  theory. 

Sensor  fusion  is  analogous  to  the  cognitive  process  used  by  humans  and  other  life  forms  to 
combine  data  from  their  senses  and  a  priori  knowledge  to  make  inferences  about  the  external  world.  One 
description  of  the  hierarchy  of  sensor  fusion  inference  ranges  from  lowest  to  the  highest:  existence  of  an 
entity,  position/velocity  determination,  identity  of  signal  sources,  behavior  of  entities,  situation 
assessment,  and  finally  threat  analysis  as  the  highest  level  of  inference[Hal92],  Thus,  this  paper  addresses 
a  relatively  low-level  sensor  fusion  inference  process. 

1.1  Background 

A  typical  networked  radar  system  consists  of  a  collection  of  remotely  located  radars.  Each  radar 
provides  a  sequence  of  time-ordered  measurements  of  position  (r,0)  in  polar  coordinate  system  or  (r,0,<|>) 
in  spherical  coordinates;  in  either  case,  the  origin  of  the  measurement  coordinate  system  is  located  at  the 
radar  antenna.  The  basic  detection  processing,  signal  processing,  and  post-detection  processing  all  occur 
at  the  radar  site.  The  digitized  position  reports,  often  called  plots,  are  forwarded  to  a  central  track 
processor  over  a  digital  data  link.  Accurate  data  registration  is  a  prerequisite  for  satisfactory  track 
initiation  and  plot-to-track  correlation;  improved  registration  reduces  requirements  for  human-machine 
interfacing  required  to  resolve  initiation  and  correlation  errors. 

At  the  central  tracking  processor,  plots  from  the  multiple  radars  update  existing  system  tracks  or 
initiate  new  system  tracks  as  appropriate.  Specifically,  the  central  tracking  processor  must  perform  the 
following  functions: 
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1 .  Transformation  of  the  radar  plots  from  local  radar  coordinates  to  system,  or  global, 
coordinates. 

2.  Correlation  or  association  of  the  radar  plots  with  the  appropriate  system  tracks.  Because 
there  are  multiple  radars  in  the  system,  more  than  one  plot  may  correlate  with  the  same  track 
over  some  arbitrary  time  interval. 

3.  Initiation  of  new  tracks  with  the  uncorrelated  plots  and  rejection  of  clutter  plots. 

4.  Track  filtering  (or  updating  with  correlated  plots)  and  track  prediction. 

5 .  Track  monitoring  and  system  track  management  (including  association  with  tracks  from 
external  sources). 

Functions  2  and  4  represent  the  heart  of  the  traditional  data  association  and  tracking  problem.  However, 
before  either  of  these  processes  can  occur  successfully,  function  1  must  be  performed;  that  is,  the 
individual  radar  data  must  have  a  common  coordinate  system  in  which  the  errors  due  to  site  uncertainties, 
antenna  orientation,  and  improper  calibration  of  range  and  time  are  minimized  so  they  do  not  cause  a 
significant  degradation  of  the  system  operation.  The  process  of  ensuring  the  required  “error  free”  (or, 
more  precisely,  controlled  error)  coordinate  conversion  of  radar  data  is  called  registration.  Thus, 
registration  is  an  absolute  prerequisite  for  multiple  radar  tracking  or  sensor  networking  in  general[Bar90]. 

1.2  Research  Motivation 

Electronic  position  location  systems  such  as  the  global  positioning  system  (GPS)  or  commercial 
satellite  survey  systems  can  locate  a  position  on  the  earth’s  surface  to  within  a  maximum  error  of  about  30 
meters  [Bar90],  This  accuracy  seems  adequate  for  radar  systems  where  the  standard  deviation  of  the  range 
measurement  error  is  greater  than  100  m.  However,  in  the  event  of  a  degraded  or  absent  GPS  signal,  a 
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mobile  ground  radar  may  have  a  positional  error  on  the  order  of  the  range  measurement  error.  This  fact 
motivates  our  analysis  of  the  specific  effect  of  radar  positional  uncertainties  (biases)  on  tracking  errors. 

The  other  systematic  error  sources  addressed,  the  range  and  azimuth  offsets  (or  biases),  are 
determined  in  part  by  radar  detection  technique  (Doppler,  mono-pulse,  etc.),  oscillator  (clock)  stability 
and  jitter,  wind  torque  and  mechanical  imperfections  (for  a  physically  rotated  antenna),  and  calibration  or 
alignment  imperfections,  etc..  By  examining  the  effect  of  the  net  position  and  range  and  azimuth  offsets, 
a  budget  may  be  obtained  that  determines  the  allowable  tolerance  of  each  of  the  above  error  sources. 

1.3  Scope  of  Research 

Of  the  sources  of  registration  error,  three  sources  which  have  proven  to  be  major  problems  in  air 
defense  and  air  traffic  control  systems  are  analyzed:  (1)  position  of  the  radar  with  respect  to  the  system 
coordinate  origin;  (2)  alignment  of  the  antennas  with  respect  to  a  common  North  reference  (that  is,  the 
azimuth  offset);  and  (3)  range  offset  errors.  Other  errors  may  exist  in  current  radar  systems;  however, 
they  have  not  been  significant  problems  in  the  past[Bar90], 

The  fourth  source  of  error,  the  inherent  inability  of  2D  radars  to  produce  the  correct  ground  range 
for  conversion  to  Cartesian  coordinates,  is  not  considered.  This  error  is  not  random,  it  always  results  in 
an  overestimate  of  the  ground  range;  the  magnitude  of  the  error  depends  on  the  aircraft  range  and 
elevation  angle.  The  only  real  solution  to  this  problem  is  the  use  of  3D  radars.  Otherwise,  the  best  that 
can  be  done  is  to  include  this  error  as  a  component  of  the  range  measurement  error.  To  remove  this 
source  of  error  from  our  model,  the  target  path  is  generated  in  three  dimensions  and  the  required 
translation  and  rotations  are  performed  by  approximating  the  Earth’s  geoid  shape  with  an  ellipsoid.  See 
Appendix  2  for  the  development  of  these  spatial  transformations.  The  range  and  azimuth  measurements 
are  then  extracted,  effectively  removing  the  ground  range  error. 
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1.4  Document  Overview 


Chapter  two  addresses  pertinent  requirements  for  determining  adequate  registration.  First  an 
analytic  framework  is  developed  to  analyze  this  issue,  with  the  Kalman  filter  model  as  the  underlying 
structure.  Background  theory  for  the  Kalman  filter  is  included  in  Appendices  one  and  two.  With  this 
analysis  model,  specific  quantitative  requirements  for  registration  are  then  developed.  Based  on  a  Chi- 
square  error  distribution,  bounds  for  the  position,  range,  and  azimuth  errors  are  established. 

Chapter  three  describes  the  development  of  the  simulation  model  and  provides  block  diagrams  for 
each  module.  In  addition,  the  operation  of  each  module  is  described  in  detail  and  the  associated 
assumptions  are  listed. 

Chapter  four  provides  the  results  of  the  simulations  and  the  accompanying  analysis.  The  target 
paths  and  representative  error  ellipse  geometries  are  included  to  help  visualize  the  results. 

Chapter  five  provides  a  summary  of  this  research  effort  and  offers  suggestions  for  further 
development. 
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2.  Summary  of  Current  Knowledge 


To  combine  information  from  multiple  sensors,  radars  in  our  case,  a  tracking  algorithm  is  often 
implemented  to  provide  measurement  accuracies  (or  statistics).  The  word  tracking  implies  that  there 
exists  a  state  estimate  S*  together  with  a  covariance  matrix  P*  for  each  detected  aircraft.  The  pair  (S*, 
P*)  could  be  obtained,  for  example,  from  the  standard  Kalman  filter  for  a  constant  velocity  plant  model 
(that  is,  no  acceleration).  The  Kalman  filter  used  in  this  application  is  a  four-state  discrete-time  filter.  If 
not  familiar  with  Kalman  filters,  refer  to  Appendices  1  and  2  for  background  information  on  Kalman 
filtering  described  in  the  context  of  tracking.  While  this  paper  is  not  intended  to  be  a  treatise  on  Kalman 
filtering,  and  definitely  doesn’t  qualify  as  one,  as  will  be  seen,  the  Kalman  filter  provides  the  tracking 
error  statistics  necessary  for  the  analysis  of  registration  errors. 

To  associate  observations  with  existing  tracks,  a  track  updating  process  typically  begins  with  a 
gating  procedure  that  is  used  to  eliminate  unlikely  observation-to-track  pairings.  Processing  is  done  at 
each  scan  using  only  data  received  on  that  scan  to  update  the  results  of  previous  processing.  This  process 
assigns  observations  to  existing  tracks  in  a  manner  that  minimizes  some  overall  distance  criterion.  At 

scan  k-1,  the  filter  forms  the  prediction  S^,  of  the  state  vector  for  use  at  time  kT.  The  measurement  at 

scan  k  is  Zk  =  HSk  +  Vk,  where  H  is  the  observation  matrix  and  V  is  zero-mean,  white  Gaussian 
measurement  noise  with  covariance  matrix  R.  The  vector  difference  between  the  measured  and  predicted 
quantities,  v  k  =  Zk  -  H  •  S^., ,  is  defined  as  the  residual,  or  innovation,  with  residual  covariance  matrix, 

®=H  P  Ht  +  R ,  where P  is  the  one-step  prediction  covariance  matrix.  The  time  subscripts  will  now 
be  dropped  for  notational  convenience.  Assume  that  the  measurement  is  of  dimension  M.  Then  defining 

d2  to  be  the  norm  of  the  residual  (or  innovation)  vector,  d2  =  v  T  0“'  -v  ,  the  M-dimensional  Gaussian 
probability  density  function  for  the  residual  is 
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e 


(2.1) 


f (v)  = 


2 


where  |©|  =  the  determinant  of  0[Bla86], 


In  either  of  the  special  cases  where  the  probability  of  detection  is  unity  or  there  are  no  expected 
extraneous  returns,  the  gate  size  should  be  infinite  for  optimal  correlation  performance.  However,  since 
one  of  the  primary  purposes  of  gating  logic  is  to  reduce  the  number  of  observation-to-track  pairings  that 
must  be  considered,  a  finite  gate  size  would  be  appropriate  even  in  these  cases.  Also,  for  non-unity 
probability  of  detection  or  during  the  presence  of  extraneous  returns,  an  optimal,  non-infinite  gate  size 
can  be  defined. 


Since  d2  is  the  sum  of  squares  of  M  independent  Gaussian  random  variables,  it  has  by  definition  a 
Chi-square  (yj)  distribution,  thus  a  correlation  gate  (G)  can  be  defined  such  that  the  plot-to-track 
association  or  correlation  decision  is  based  on  a  Chi-squared  test  of  the  following  form: 

Pp  ~z]r[pp  +PZ]"' [sp  -z]<  G  (2‘2) 


where  Sp  denotes  the  position  components  of  S*  extrapolated  to  the  time  at  which  the  next  measurement 
Z  is  obtained;  that  is, 


S  =  <E(  At)S  *  P  =<t(A?)P*<t(A?)r 


(2.3) 


where  <J(A/)  is  the  state  transition  matrix  for  a  time  At.  Also  Z  is  the  measured  position  (in  2D  Cartesian 
coordinates),  Pp  the  covariance  submatrix  of  P  for  the  position  components,  and  Pz  the  covariance 
matrix  for  the  measurement. 
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If  a  radar  measurement  Z  from  radar  A  satisfies  the  gate  test  defined  by  (2.2),  then  Z  is  used  to 
update  the  track  through  the  estimation  procedure.  If  more  than  one  measurement  from  radar  A  satisfies 
the  gate  test,  then  an  ambiguity  resolution  logic  is  necessary  to  select  one  measurement  for  the  track 
update  process.  This  can  be  accomplished  with  an  optimal  assignment  algorithm  such  as  the  Munkres 
algorithm.  However,  if  no  measurement  from  radar  A  satisfies  the  gate  test,  then  the  gate  G  may  be 
enlarged  by  adding  a  “maneuver  term”: 

G'=G +  (\/C)(AM)2  ,AM  =  At2 12  (2.4) 

for  a  maneuver  or  acceleration  factor  A.  The  normalizing  factor  C  can  be  defined  as  either  the  minimum 
eigenvalue  of  the  joint  covariance  matrix  [Pp  +  Pz]  or  the  nth  root  of  the  determinant  of  the  (n  x  n) 
covariance  matrix.  If  the  correlated  plot  is  in  the  maneuver  gate  (but  not  in  the  nonmaneuver  gate)  for 
two  or  more  successive  scans,  then  a  “maneuver”  might  be  declared  and  a  special  set  of  maneuver  gains 
or  filter  constants  used  to  update  the  track.  A  maneuver  gate  is  not  implemented  in  the  model  since  only 
linear  target  paths  are  generated.  This  is  done  to  simplify  the  Kalman  filter  algorithm.  The  Kalman  filter 
is  the  optimal  estimator  for  a  linear  path,  however  a  maneuvering  target  requires  an  extended  Kalman 
filter  in  which  its  state  vector  is  augmented  with  additional  components. 

Referring  to  Figure  2-1,  suppose  that  an  aircraft  is  tracked  by  two  radars  denoted  as  radars  A  and 
B  whose  detectable  coverage  areas  overlap.  Where,  or  even  if,  a  plot  from  radar  B  falls  in  the  correlation 
gates  depends  both  on  the  random  measurement  errors  of  the  radar  and  the  magnitude  of  the  position, 
range,  and  azimuth  biases  between  the  two  radars. 
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Radar  B 


Figure  2-1  Constant  Measurement  Error  Ellipses 

If  there  are  no  registration  errors  or  biases,  then  the  plot  from  radar  B  should  fall  within  the 
nonmaneuver  gate  G  most  of  the  time.  Presumably,  gate  G  was  chosen  to  ensure  that  plots  from  a 
nonmaneuvering  aircraft  will  satisfy  (2.2)  with  a  probability  in  the  range  0.90  to  0.99,  based  on  the 
characteristics  of  the  random  errors  and  the  actual  system  design  goals.  Similarly,  if  the  biases  are  small 
with  respect  to  the  random  errors,  then  the  plot  from  radar  B  should  be  in  the  nonmaneuver  gate  most  of 
the  time,  although  the  exact  probability  will  be  less  than  the  design  goal.  Similarly,  if  the  biases  are  small 
with  respect  to  the  random  errors,  then  the  plot  from  radar  B  should  be  in  the  nonmaneuver  gate  most  of 
the  time,  although  the  exact  probability  will  be  less  than  the  design  goal. 

On  the  other  hand,  if  the  biases  are  relatively  large  with  respect  to  the  random  errors,  perhaps 
approximately  of  the  size  of  the  gate  G  or  even  G',  then  the  plot  may  fall  between  the  nonmaneuver  and 
the  maneuver  gates: 


G<[sp-Zf[pp+Pz]_,[sp-Z]<G' 


(2.5) 
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Although  it  is  unlikely  that  a  tracking  system  would  be  designed  to  declare  a  maneuver  on  one  maneuver 
gate  correlation,  the  possibility  now  exists  that  a  maneuver  could  be  falsely  declared  if  radar  A 
subsequently  fails  to  detect  the  aircraft  on  the  next  scan.  If  the  plot  from  radar  B  is  used  to  update  the 
track,  then  the  bias  is  superimposed  on  the  state  estimate,  with  a  loss  of  system  track  accuracy.  If  the  plot 
is  simply  discarded,  then  the  system  may  have  a  delayed  response  to  an  actual  aircraft  maneuver; 
certainly,  there  is  a  loss  of  information.  Finally,  if  the  biases  are  very  large  with  respect  to  the  random 
errors,  then  the  plot  from  radar  B  will  not  correlate  with  the  track  at  all,  in  which  case  the  system 
eventually  will  initiate  a  second  track  for  the  same  aircraft. 

2.1  Random  versus  Systematic  Errors 

It  was  asserted  in  the  preceding  paragraph  that  the  biases  or  registration  errors  in  plot  data  would  degrade 
system  track  accuracy  if  used  to  update  a  system  track.  However,  this  is  only  one  example  of  a  more 
general  problem  in  the  general  theory  of  estimation.  To  consider  the  more  general  problem,  let  X  be  a 
random  variable  with  an  expected  value  p  and  standard  deviation  <7;  that  is, 

ft  =  E[X],  ct2  =  E[(X  -  p)2]  (2.6) 

A  random  sample  x  from  X  can  be  represented  as 

x  =  p  +  e  (2.7) 

where  s  is  the  zero-mean  random  component  of  X;  therefore, 

E[e]  =  0,  ct2  =  E[e2]  (2.8) 

In  many  radar  tracker  designs,  it  is  assumed  that  the  radar  measurements  have  the  properties  just 
outlined.  In  particular,  it  is  assumed  that  p  represents  the  true  value  of  the  aircraft  position.  Tracking  or 
filtering  is  a  process  for  estimating  the  true  aircraft  position  p  from  a  sequence  of  measurements;  that  is, 
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random  samples,  taken  over  time.  Thus,  if  the  set  of  measurements  {x,,  x2,  x3,.„,  xN}  is  available;  then 
the  estimate  x*  of  p  is  given  by 

X*  =  ajxi  +  a2x2  +  a3x3  +. . .  +  aNxN  (2.9) 

where  the  set  of  nonnegative  scalars  {a],  a2,  a3, .  .  . ,  aN},  determined  by  the  gains  of  the  filter  algorithm, 
satisfies  the  usual  convexity  criterion: 


a1  +  a2  +  a3  +  ...  +  aN=l  (2.10) 

If  it  is  true  that  E[e]  =  0,  then  x*  is  an  unbiased  estimate  of  the  true  position  p  of  the  aircraft. 
However,  if  E[s]  =  p,  where  p  is  not  zero  (P  *  0),  then  by  (2.10)  it  follows  that 

E[x*]  =  p  +  p  (2.11) 

and  the  estimate  x*  is  a  biased  estimate  of  p.  Moreover,  the  mean  square  error  (mse),  V,  of  a  biased 
estimate  of  x*  is  larger  than  an  unbiased  estimate;  that  is, 

V  [x*]  =  E  [(x*  -  p)2]  =  (a,2  +  a22  +  a32  +  . . .  +  aN2)a2  (2.12) 

if  E[e]  =  0;  whereas 

V[x*]  =  E[(x*  -  p)2]  =  (a,2  +  a22  +  a32  +  . . .  +  aN2)cr2  +  p2  (2.13) 

if  E[e]  =  p. 

From  this  it  follows  that  the  process  of  filtering  or  state  estimation  “averages”  the  random 
measurement  errors  to  reduce  the  variance  or  mse  of  the  estimate  of  the  true  state  p.  However,  the 
filtering  process  cannot  remove  or  even  reduce  the  magnitude  of  the  bias  or  systematic  errors. 
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In  the  context  of  a  multiple  radar  tracking  system,  the  presence  of  registration  errors  will  result  in 
a  track  mse  larger  than  that  which  should  be  achievable  theoretically.  If  the  registration  errors  are 
sufficiently  large,  then  multiple  radar  tracking  will  be  less  accurate  than  single  radar  tracking.  In  the 
worst  case,  registration  errors  can  result  in  a  failure  to  correlate  multiple  radar  measurements  from  a 
common  aircraft  with  a  common  track,  either  reducing  the  system  effectively  to  a  single  radar  system  or 
even  leading  to  multiple  tracks  for  the  same  aircraft.  Because  the  fundamental  objective  of  multiple  radar 
tracking  is  to  ensure  track  continuity  for  an  aircraft  as  it  moves  through  multiple  radar  coverage 
envelopes,  registration  errors  can  defeat  the  very  purpose  of  a  multiple  radar  system.  Therefore,  the  basic 
requirement  for  registration  is  to  ensure  that  plots,  that  is,  radar  measurements,  from  a  common  aircraft 
will  be  in  the  nonmaneuver  correlation  gate  (in  the  absence  of  maneuvers)[Bar90]. 

2.2  Quantitative  Requirements  for  Registration 

At  this  point,  the  need  for  radar  registration  should  be  obvious.  The  next  question  therefore  is, 
how  well  must  radars  be  registered?  Before  it  will  be  possible  to  address  this  question  directly,  some 
results  from  the  distribution  theory  for  normally  distributed  random  variables  must  be  stated.  Based  on 
this  theory,  an  analysis  model  for  the  effects  of  registration  errors  on  plot  correlation  will  be  developed. 
Finally,  this  model  will  be  applied  to  derive  some  quantitative  requirements  for  the  three  major  sources  of 
registration  error:  sensor  position,  range  offset,  and  azimuth  offset.  These  results  are  discussed  at  length 
in  many  textbooks  on  multivariate  statistical  analysis. 

For  this  discussion,  assume  that  {Xk  |  k  =  1,  2, .  . . ,  N}  is  a  set  of  normally  distributed,  scalar 
random  variables  with 


E[Xk]  =  pk,  E[(Xk  -  pk)2]  =  ok2 
that  is,  each  Xk  is  distributed  as  a  N(pk,  ak)  random  variable. 

Now  consider  a  random  variable  Z  defined  as  the  normalized  sum  of  squares  of  the  Xk: 


(2.14) 
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z=Z 

*=n 


Xy 


Mk 


^2 


V 


(2.15) 


The  random  variable  Z  is  distributed  as  a  Chi-squared  random  variable  with  N  degrees  of  freedom, 
denoted  by  x2(N).  The  distribution  x2(N)  is  often  called  the  central  Chi-squared  distribution  because  each 
of  the  normalized  terms, 

(2.16) 

°k 


from  the  sum  of  (2.15)  is  distributed  as  aN(0,  1)  random  variable,  which  is  often  noted  X'k~N(0,l). 
Equivalently,  each  normalized  term  squared 


z»=(x;)!  = 


V  °k  J 


(2.17) 


is  a  x2(l)  random  variable.  If  the  means  pk  are  omitted  from  the  sum  in  (2.15),  that  is. 


N 

Z=  I 

k=\' 


X y 


\Uk  J 


(2.18) 


then  Z  is  distributed  as  a  noncentral  Chi-squared  random  variable  with  parameter  X,  denoted  by  Z  ~  x2(N, 
X).  The  noncentrality  parameter  X  is  defined 


N 

k=\' 


ru 
Mk 

\akJ 


(2.19) 


The  general  theory  just  outlined  can  be  extended  to  random  vectors  in  a  straightforward  manner. 
For  this  case,  let  X  denote  a  normally  distributed  random  vector  in  RN  with  mean  p  and  covariance  matrix 
P.  Then  the  quadratic  form  4  defined  by 
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^  =  (X  -  |x)TP''(X  -  p.) 


(2.20) 


is  distributed  as  a  x2(N)  random  variable,  q  ~  x2(N) .  Similarly,  if  the  quadratic  form  £,  is  defined 

4  =  XTP-'X  (2.21) 

then,  %  ~  x2(N,  X)  ,  where  the  noncentrality  parameter  X  is  given  by 

X  =  p.TP‘V  (2.22) 


2.3  Analysis  Model 

With  the  help  of  the  general  theory  just  outlined,  it  is  now  possible  to  develop  a  mathematical 
model  with  which  the  quantitative  effects  of  registration  errors  can  be  examined,  specifically,  the  effects 
on  multiple  radar  tracking  and  correlation.  To  proceed,  let  Z  denote  a  measurement  of  an  aircraft  position 
from  a  radar  in  the  system.  We  assume  that  tracking  is  performed  in  a  Cartesian  coordinate  system,  either 
R4  or  R6,  (the  positions  x,y  or  x,y,z  and  their  corresponding  velocities).  Consequently,  it  may  be  assumed 
that  Z  is  in  system  coordinates;  that  is,  the  measurement  vector  Z  is  of  the  form  [x,  y]  or  [x,  y,  z],  rather 
than  the  natural  radar  polar  coordinates  [r,  0]  or  [r,  0,  cp]. 

Although  it  cannot  be  proven  rigorously,  we  may  assume  that  the  measurement  Z  is  a  normally 
distributed  random  vector  with  the  mean  equal  to  the  true  aircraft  position.  This  assumption  is  partially 
justified  for  two  reasons.  First,  the  radar  range  measurement  r  generally  follows  a  Rayleigh  distribution. 
The  azimuth  measurement  0  certainly  is  not  uniformly  distributed  over  the  interval  [0,  2 n);  however  the 
value  of  the  azimuth,  in  effect,  is  quantized  in  the  radar  signal  processor  and  in  the  analog-to-digital 
conversion  process.  Consequently,  it  is  reasonable  to  assume  that  0  is  approximately  uniformly 
distributed  over  some  subinterval  of  [0,  27t)[Bar90].  It  can  be  shown  that  the  polar-to-Cartesian 
coordinate  conversion 
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x  =  r  cos(0),  y  =  r  sin(0) 


(2.23) 


yields  two  independent  N(0,  1)  random  variables  if  r  follows  a  Rayleigh  distribution  and  0  is  uniformly 
distributed  over  [0,  2n)[Pap91], 

If  the  random  variables  x  and  y  are  ~N(0,c)  and  independent. 


fxr(x>y) 


1  e~(x1+y1)l2<r1 

2  na1 


then  to  find  the  probability  distribution  function  in  polar  coordinates, 


(2.24) 


fR®(r,0)  = 


fxA^y) 

\Ax,y)  | 


(2.25) 


where  |J(x,y)|  is  the  absolute  value  of  the  Jacobian  of  the  polar  to  Cartesian  transformation  which  is, 


J{x,y) 
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where  |..|  is  the  determinant  of  the  matrix.  From  (2.24),  it  can  be  concluded  that, 


fR@(r,0)  =  r-  fXY{x,y)  = 


2n<j 2 


,r> 


0,|#|  <  n 


(2.26) 


(2.27) 


and  0  otherwise.  This  is  a  product  of  a  function  of  r  times  a  function  of  0.  Hence  the  random  variables  r 
and  0  are  independent  with 


fR(r)  =  -^1e~rl'la\fQ{6)  =  - — ,r  >  0,|^|  < 

<7  2 7Z 


(2.28) 


15 


Where  the  proportionality  constant  is  chosen  to  make  the  area  (total  probability)  of  each  term  equal  to  1. 
Thus  it  follows  that  if  x  and  y  are  ~N(0,g)  and  independent,  the  random  variables  r  and  0  are  independent 
with  r  being  Rayleigh  distributed  and  0  being  uniformly  distributed. 

Secondly,  empirical  measurements  gathered  by  Hughes  Aircraft  Company  over  past  years 
suggests  that  the  random  errors  in  Z  are  at  least  approximately  normally  distributed[Bar90],  It  is  true, 
however,  that  there  is  a  distinct  correlation  in  the  errors  between  x  and  y.  This  is  the  reason  for 
employing  a  Kalman  filter  capable  of  tracking  in  two  dimensions.  This  allows  us  to  take  advantage  of  the 
coupling,  or  correlation,  of  the  x  and  y  measurements.  Recall  that  the  actual  radar  measurements  are  in 
polar  coordinates,  and  any  change  in  range  or  azimuth  generally  results  in  a  corresponding  change  in  x 
and  y.  Of  course,  in  the  special  case  of  a  target  path  co-linear  with  either  the  x  or  y  axis,  a  change  of 
range  will  not  result  in  changes  in  both  x  and  y. 

Assume  that  a  time-ordered  sequence  of  correlated  measurements  {Zh  Z2,  Z3, .  . . ,  Zk)  are 
processed  by  a  Kalman  filter  to  obtain  a  track,  which  consists  of  a  state  estimate  S*  and  a  covariance 
matrix  P*  for  the  state  estimate.  Therefore,  the  state  estimate  S*  is  a  linear  combination  of  the 
measurements  {Z1;  Z2,  Z3, . . . ,  Zk}.  With  the  assumption  of  normality  for  the  errors  in  Z,  and  the  natural 
assumption  of  independence  of  the  time  sequence  of  measurements,  it  follows  that  the  state  estimate  S*  is 
a  normally  distributed  random  vector. 

Now  consider  the  correlation  criterion  given  in  (2.2) 


(2.29) 


The  quadratic  form  £,  in  (2.24)  is  distributed  as  a  x2(M,  X)  random  variable  with  the  number  of  degrees  of 
freedom  M  equal  to  the  dimension  of  the  measurement  vector  Z. 
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In  theory,  the  non-centrality  parameter  X  should  be  zero  if  the  measurement  Zk+]  was  obtained 
from  the  same  aircraft  represented  by  the  track  vector  S*.  The  parameter  X  can  be  other  than  zero  if 
either  (1)  the  measurement  Zk+1  was  obtained  from  a  different  aircraft  than  that  represented  by  the  track  or 
(2)  there  are  biases  that  would  create  an  apparent  difference  without  the  random  measurement  errors.  For 
this  application,  X  will  represent  the  total  normalized  biases  in  the  measurement  vector  Z;  that  is, 

X=br[F,+PzI'b  <2-M> 

Note  that  (2.25)  assumes  that  the  measurement  vector  Z  is  of  the  form 

Z=p+b+s  (2.31) 

where  p  is  the  true  aircraft  position  and  8  is  the  vector  of  random  errors.  Assuming  that  the  mean  of  the 
prediction  Sp  is  p,  then 


E\ Sp-Z]=b  (2.32) 

Equations  (2.29)  and  (2.30)  are  the  model  with  which  the  impact  of  registration  errors  can  be 
analyzed  quantitatively.  The  non-maneuver  gate  G  and  the  maneuver  gate  G'  are  chosen  to  obtain  a 
specified  probability  of  correlation  of  plots  that  represent  the  same  aircraft  as  the  track.  For  example,  G  is 
chosen  from  a  y2(M)  probability  distribution  to  satisfy 

Prob[£,  <  G]  >  p0  (2.33) 

The  objective  now  is  to  define  an  error  “budget”  for  the  sources  of  registration  error  such  that 

Prob[^  <  G]  >  p0  -  Ap  (2.34) 
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where  the  correlation  statistic  4  is  distributed  as  a  x2(M,  X)  random  variable  with  parameter  X  given  by 
(2.30),  and  Ap  >  0  is  the  reduction  in  the  correlation  probability  that  can  be  tolerated  if  the  system  still  is 
to  meet  the  system-level  requirements  for  tracking  or  track  accuracy. 

2.4  Quantitative  Requirements 

To  develop  some  specific  requirements  for  registration,  consider  first  the  correlation  gate  size  G. 
One  common  rule  of  thumb  in  tracking  systems  is  to  choose  G  such  that 

Probg  <G]  =  0.99  (2.35) 

If  the  measurement  vector  Z  is  in  R2,  then  the  inverse  Chi-square  with  N  =  2  degrees  of  freedom  and  a 
probability  of  .99  yields  G  =  9.2;  similarly  for  Z  in  R3,  G  =  1 1 .3.  A  correlation  probability  of  0.99  may 
be  an  excessive  requirement  considering  the  probability  of  detection  of  many  surveillance  radars,  which 
often  are  specified  to  be  only  0.8  to  0.9[Bar90],  Consequently,  a  correlation  probability  of  0.95  would 
seem  adequate  for  most  applications.  Given  that  0.95  is  adequate,  then  from  Figure  2-2,  a  gate  size  of  9.2 
allows  a  total  bias  parameter  X=\2  for  the  two-dimensional  case  or  X  =  1 .5  for  the  three-dimensional 
case. 


Figure  2-2  Non-central  Chi-square  (x2)  cdf  vs.  non-centrality  parameter 
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The  three  major  sources  of  registration  error  are  sensor  position,  range  offset,  and  azimuth  offset. 
The  next  step  in  the  analysis,  therefore,  is  to  find  the  maximum  error  in  each  variable  consistent  with  the 
bound  of  the  parameter  X.  In  the  following  discussion,  only  the  case  of  a  two-dimensional  measurement 
vector  will  be  considered  since  the  model’s  measurements  are  in  those  dimensions. 

Consider  an  aircraft  tracked  by  radars  A  and  B.  For  convenience,  let  PA  and  PB  denote  the 
covariance  matrices  for  the  measurement  errors  at  radars  A  and  B,  respectively.  The  registration  problem 
is  to  ensure  that  the  bias  vector  b  is  sufficiently  small  that  the  measurement  from  radar  B  will  correlate  in 
the  non-maneuver  gate  with  the  system  track. 

That  is,  the  problem  now  is  to  bound  the  errors  that  contribute  to  the  bias  vector  b  from  the 
inequality 


Pp+P 


<1.2 


(2.36) 


Obviously,  it  would  be  convenient  to  eliminate  one  of  the  variance  parameters.  Many  tracking  systems, 
in  an  attempt  to  maintain  sensitivity  to  possible  aircraft  maneuvers,  bound  the  gains  in  the  Kalman  filter 
from  below[Bar90],  One  rule  of  thumb  is  to  allow  the  gains  to  decrease  to  the  level  that  produces  steady 
state  position  variances  to  be  approximately  50  percent  of  the  corresponding  measurement  variances. 
Thus,  (2.36)  becomes 


br  [0.5  •  P^  +  PB  ]_1  b  <  1 .2 


(2.37) 


Before  continuing  further,  consider  briefly  the  nature  of  the  radar  measurement  error.  In  the  radar 
measurement  plane,  that  is,  the  [r,  0]  plane,  the  covariance  matrix  P  is  defined 


P'  - 

rR 


V 


2  2 


r  <j, 


ej 


(2.38) 
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where  crr  and  ct0  are  the  standard  deviations  of  the  range  and  azimuth  measurements,  respectively;  and 
where  r  is  the  range  from  the  radar  to  the  aircraft.  In  modem  surveillance  radars,  the  range  measurement 
often  is  much  more  accurate  than  the  cross-range  measurement,  which  can  result  in  highly  elliptical 
equal-probability  contours  (or  error  ellipses)  in  the  plane.  The  error  ellipses  in  R2  have  major  and  minor 
semi-axes  equal  to  the  square  roots  of  the  two  diagonal  components  of  the  matrix  in  equation  (2.38). 

The  basic  concept  of  an  error  ellipse  describes  the  eccentricity  of  a  single  radar  error  contour 
determined  by  the  ratio  of  the  range  error  ar  to  the  cross-range  error  r-ae.  In  this  example,  the  range  error 
a,  is  smaller  than  the  cross-range  error  r-ae,  which  creates  highly  elliptical  errors  in  the  plane  for  both 
radars.  In  Figure  2-1,  the  sum  is  nearly  circular  because  of  the  orthogonal  geometry  of  the  aircraft  and  the 
two  radars. 

The  conversion  from  radar  coordinates  to  the  Cartesian  plane  rotates  an  error  ellipse  in  the  plane 
by  a  unitary  matrix  of  the  form 


^cos#  -sin#N 
vsin#  cos#  J 


(2.39) 


where  9  is  the  angle  to  the  target,  measured  counterclockwise  from  the  abscissa,  (the  x  axis),  of  the  radar 
coordinate  system.  Thus, 


Pfl=U-P'-Ur  (2-40) 

for  radar  R.  Note  that  the  unitary  transformation  U  merely  rotates  an  error  ellipse  with  respect  to  the  axes 
of  the  Cartesian  plane;  it  does  not  change  the  basic  shape  or  area  of  an  error  ellipse.  Because  the  angle  0 
will  be  different  for  different  radars  viewing  the  same  aircraft,  the  orientation  of  the  error  ellipses  PA  and 
PB  can  be  quite  different,  even  if  the  error  matrices  P'A  and  P'B  in  the  respective  radar  planes  are  similar. 
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Now,  consider  the  problem  of  a  bias  b  due  to  an  error  in  the  knowledge  of  the  position  of  radar  B 
with  respect  to  radar  A.  Equation  (2.37)  bounds  the  magnitude  of  b  relative  to  the  sum  of  the  tracking 
and  measurement  error  variances  [0.5-PA+PB].  The  worst  case  geometry  occurs  when  (1)  the  aircraft 
under  consideration  is  at  the  midpoint  of  the  line  segment  that  joins  the  two  radar  locations  and  (2)  the 
bias  vector  b  is  parallel  to  that  line. 


b2  =  b'b 


:  1.8  •  min  cr2  (a), cr2  (5)]  =1.8-  a 2  (min) 


(2.41) 


Now  the  limit  of  the  allowable  absolute  position  bias  is  defined  as, 

|by,|  =  1.34-  crr  (min)  (2.42) 

This  assumes,  of  course,  that  the  range  error  is  no  greater  than  the  cross-range  error.  However  if  this  is 
not  the  case,  then  the  cross-range  error  for  some  minimum  detection  range  could  be  used  instead.  In  this 
case,  (2.37)  becomes, 


b2  =|brb|<1.8'min[r2  ■a2e{A),r2  •cr^(5)]=  1.8-r2(7^(min)  (2-43) 

where  Ge(A)  denotes  the  standard  deviation  of  the  azimuth  measurements  from  radar  A,  and  r  is  the 
down-range  distance  from  the  radar  to  the  target.  Therefore,  taking  the  square  root  of  equation  (2.43) 
gives  the  maximum  tolerance  bP  on  the  magnitude  of  the  position  error, 

Jbp  |  =  1.34-r -cr^(min)  (2.44) 

where  ae(min)  is  the  standard  deviation  of  the  radar  with  the  smallest  azimuth  measurement  error. 

Next,  consider  the  problem  of  range  offset  errors.  The  worst  case  geometry  is  the  same  as  that  for 
radar  position  errors.  The  range  bias  problem,  however,  is  not  a  relative  problem;  that  is,  each  radar  can 
contribute  an  error.  Whether  these  errors  are  additive  or  partially  cancel  each  other  depends  on  the  exact 
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geometry  as  well  as  the  sign  of  the  bias.  For  the  geometry  under  consideration,  the  worst  case  is  for 
errors  of  the  same  sign.  Therefore,  the  equivalent  form  of  (2.44)  is  the  following: 

|2b*|  =  1.34-ov  (min)  (2-45) 

where  bR  is  the  maximum  tolerance  for  the  range  bias  at  each  radar;  or  equivalently, 

1^1  =  0.67-07  (min)  (2-46) 

Lastly,  consider  the  problem  of  antenna  alignment  or  azimuth  offset  errors.  The  error  or  offset  b 
of  the  measurement  from  radar  B  will  be  parallel  to  the  cross-range  error  (due  to  azimuth  measurement 
errors).  The  magnitude  of  the  error,  be,  is  given  by 

b,=r*(A&)  (2.47) 


where  Ap9  is  the  azimuth  offset. 

The  geometry  in  this  case  is  more  complex;  however,  some  empirical  observations  will  suffice  to 
show  that  the  worst  case  occurs  for  an  orthogonal  aircraft-radar  geometry.  Assuming  that  the  axes  of  the 
Cartesian  coordinate  system  are  aligned  so  that  the  origin,  aircraft,  and  the  two  radars  compose  the  four 
corners  of  a  rectangle.  The  azimuth  offset  error  will  be  along  the  cross-range  error  with  respect  to  radar 
B,  which  is  parallel  to  the  range  error  for  radar  A.  Therefore,  it  follows  that 

bg  =|brb|<1.2-[o.5-crr2(^)+r2  .<jj(fl)]  (2-48) 

If  it  is  assumed  that  cr(s4)  «  r<j0(B) ,  then  it  follows  that  a  reasonable  bound  for  the  azimuth  offset 
Ape  is  given  by 
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(2.49) 


|A^|<1.1  -g6 

Here,  the  azimuth  offset  has  been  treated  as  a  relative  error  at  radar  B,  which  is  adequate  for  a  system  of 
only  two  radars. 

However,  if  the  down-range  measurement  variance  cr  is  greater  than  the  cross-range 
measurement  variance  r-Ge,  then 

b]  =  \bT b|  <  1 .2  •  [.5  •  r2  •  <7 1 (. A)  +  r 2  •  <r|  (5)]  (2-50) 

and  using  equation  (2.47)  and  the  fact  that  the  worst-case  geometry  will  now  occur  when  the  target  is 
centered  between  the  two  radars,  and  assuming  the  two  azimuth  variances  are  equal, 

\Af3e\<4U-^y2-(je=\34-cje  (2-51) 

which  is  about  22  A>  greater  than  the  maximum  allowable  bias  using  the  previous  assumption. 

In  a  system  with  three  or  more  radars,  it  would  be  necessary  to  align  each  radar  with  respect  to 
true  north.  Therefore,  the  actual  requirement  in  a  multiple  radar  system  using  the  assumption  that  the 
cross-range  measurement  variance  is  much  greater  than  the  down-range  variance  is 

\Aj3d\<0.55-ag  (2.52) 

Assuming  the  down-range  variance  is  greater  than  the  cross-range  variance, 

|A/?,  |  <0.67- <7,  (2.53) 

The  results  of  these  single-source  error  analyses  are  summarized  in  Table  1  in  the  form  of  a 
registration  error  budget.  However,  in  actuality,  these  errors  occur  simultaneously,  resulting  in  a 
cumulative  error.  If  all  three  error  sources  are  considered  together  as  additive  vectors,  then  by  the  nature 
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of  vector  addition,  the  error  budget  must  be  reduced  by  a  factor  of  V3  ;  the  result  is  shown  in  the  right- 
hand  column  of  Table  2.1. 


Table  2-1  Registration  Error  Budgets 


Error  Source 

Single-Source  Tolerance 

Multi-Source  Tolerance 

Radar  Position 

1.34-0r(min) 

0r(A)«  r-0@(B) 

O.77-0r(min) 

0r(A)«  r-0e(B) 

1.34-r-ae(min) 

0r(A)  *  r-00(B) 

O.77-r-0e(min) 

0r(A)  *  r-00(B) 

Range  offset 

0.67-ar(min) 

O.39-0r(min) 

Azimuth  offset 

0.55-00 

0r(A)«  r-00(B) 

O.32-0e 

0r(A)«  r-0e(B) 

0.67-00 

0r(A) «  r-0e(B) 

0.39-00 

0r(A)  *  r-09(B) 

Note:  CTr(min)  is  the  minimum  standard  deviation  over  all  radars  in  the  system.  The  bound  for  the 
azimuth  bias  can  be  taken  relative  to  each  site. 


2.5  Summary 

Developed  in  this  chapter,  is  the  model  by  which  the  effects  of  registration  error  in  a  networked 
radar  system  may  be  quantitatively  analyzed.  Also  introduced  is  the  concept  of  a  correlation  gate  as  the 
sum  of  squared  Gaussian  distributed  random  variables,  constituting  a  non-central  Chi-squared  random 
variable  with  degrees  of  freedom,  N  equal  to  the  number  of  dimensions  in  the  measurement  vector.  It  is 
observed  that  the  correlation  gates  most  often  take  the  form  of  “error  ellipses”  in  two-dimensional 
measurement  space.  Finally,  the  worst-case  geometries  for  the  target,  radars  and  global  coordinate  origin 
are  presented.  The  three  sources  of  error  considered,  the  position,  range,  and  the  azimuth  biases  are  to  be 
quantitatively  analyzed  individually  and  cumulatively. 
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3.  Simulation  Model  Description 


3.1  Development  of  the  Simulation  Model 

One  of  the  goals  of  this  effort  is  the  creation  of  a  user-friendly  software  model  which  may  be  used 
in  further  analysis  of  networked  radar  systems.  Since  there  exists  a  number  of  industry  accepted  and 
validated  software  models  to  simulate  various  radar  types,  resources  are  not  expended  in  duplicating  this 
capability.  In  addition,  future  efforts  will  be  made  to  import  radar  plots  generated  from  other  tools  as  the 
input  to  this  model.  A  model  is  desired  that  will  be  easily  altered  to  allow  for  the  analysis  of  differing 
sensor  measurement  characteristics,  including  azimuth  and  range  measurement  variances.  It  is  also 
desirable  to  permit  arbitrary  placement  of  sensor  positions,  and  to  change  the  position,  range,  and 
azimuth  biases  easily. 

Because  of  these  desired  characteristics,  it  is  decided  to  implement  the  model  in  the  graphical 
simulation  environment  Simulink®,  an  extension  to  the  widely  accepted  commercial  software  package 
MatLab®.  One  implication  of  using  Simulink®  is  that  not  all  MatLab®  functions  are  directly  callable  from 
its  block  diagram  simulation  environment.  To  access  certain  MatLab®  functions,  or  custom  written 
functions,  S-functions  must  be  written.  S-functions  use  a  special  calling  syntax  that  enable  you  to  interact 
with  the  MatLab®  ordinary  differential  equation  solver. 

Another  difficulty  using  Simulink®  is  that  (m  x  n)  matrices  are  decomposed  and  reconstructed 
into  (m  n  x  1)  vectors  and  passed  from  block-to-block  during  a  simulation.  To  perform  matrix 
calculations,  it  is  advantageous  to  use  the  MatLab®  Digital  Signal  Processing  Blockset®  to  ease  the 
burden  of  resizing  these  vectors  into  their  original  dimensions. 
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Despite  these  restrictions,  Simulink®  provides  an  intuitive  environment  for  model  creation,  allows 
easy  model  parameter  changes,  and  allows  Monte-Carlo  simulations  to  be  performed  from  the  MatLab® 
command  line. 

3.2  Simulation  Model  Block  Diagram  Description 

The  block  diagrams  that  constitute  the  simulation  model  are  now  considered.  Note  that  in  the 
following  illustrations,  thicker  lines  denote  vectors  and  the  thinner  lines  scalars. 


To  Workspaces 


Figure  3-1  Top-Level  Block  Diagram 


Figure  3-1  illustrates  the  top-level  diagram  of  the  simulation  model.  At  the  extreme  left  of  the 


diagram  is  the  target  path  generation  and  coordinate  conversion  block.  This  block  passes  the  target  paths 


from  the  two  radars  to  the  Kalman  filter  and  Chi-square  test  block.  The  binary  output  of  the  Chi-square 


test  is  then  passed  to  the  track  maintenance  blocks  which  perform  sequential  testing  to  determine  the  track 
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state.  The  “to  workspace”  blocks  collect  test  data  which  is  then  accessed  from  the  MatLab®  command 
line  to  implement  Monte  Carlo  analysis. 


Figure  3-2  Path  Generation  Block  Diagram 

Figure  3-2  illustrates  the  block  diagram  responsible  for  generating  the  target  path,  measurement 
variances,  position  bias,  range  bias,  and  azimuth  bias.  The  coordinate  transformation  block  performs  a 
global  to  local  transformation,  injects  the  biases,  and  then  performs  the  local  to  global  transformation  for 
each  radar.  The  relevant  equations  for  the  coordinate  transformation  are  included  in  Appendix  three. 


Figure  3-3  Coordinate  Conversion  and  Bias  Insertion  Block  Diagram 
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Figure  3-3  illustrates  the  internal  block  diagram  responsible  for  local -to-global  coordinate 
transformation,  and  the  insertion  of  position,  range  and  azimuth  biases.  Notice  that  the  position  biases  are 
inserted  in  cartesian  coordinate  space,  while  the  range  and  azimuth  biases  are  inserted  in  spherical 
coordinate  space.  It  may  be  intuitively  obvious  that  the  nature  of  these  bias  will  influence  the  gate 
correlation  process  differently.  More  will  be  discussed  about  this  topic  in  the  analysis  of  results  section  of 
this  paper. 


Sk/k 

Sk+l/k 

Sk/k-1 

Innov. 

Pk/k 

ZJn  Pk+l/k 
Kalm.  Gain 

HPk/k-lH* 
H*Sk/k-l 
Hk*  Sk+l/k 


Filter  Block 


Z  in 


Sk/k 
Sk+l/k 
Sk/k-1 
Innov. 

Pk/k 
Pk+l/k 
Kalm.  Gain 
Rb 

Pk/k- 1 
H*Sk/k-l 
Hk*Sk+l/k 


Filter  Block 


dA2 

+ 

bias 

Chi2_test 

I _ 

Gate  Test 


dA2 


Chi2  Test 


Figure  3-4  Kalman  Filter  and  Chi-square  Gate  Test  Block  Diagram 


Figure  3-4  illustrates  the  block  diagram  responsible  for  generating  the  target  tracks  and 
accompanying  statistics  which  are  passed  to  the  Chi-square  gate  test  blocks.  Notice  that  the  bias  is  simply 
the  difference  of  the  state  position  components  from  radar  one  and  the  measurements  from  radar  two. 

The  Kalman  filters  also  pass  the  required  measurement  error  covariance  and  state  prediction  error 
covariance  required  by  the  Chi-square  gate  test  blocks. 
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Figure  3-5  Kalman  Filter  State  Estimation  Loop 


The  Kalman  filter  consists  of  two  distinct  recursive  “loops”.  Figure  3-5  represents  the  Kalman 
filter  state  estimation  loop  which  recursively  estimates  the  current  state  and  also  calculates  the  one-step 
prediction  of  the  next  measurement.  It  can  be  seen  that  the  prediction  is  delayed  by  the  period  between 
individual  measurements  which  corresponds  to  the  scan  time  for  the  mechanically  rotated  radar.  The 
measurement  components  of  the  predicted  state  vector  are  extracted  and  then  subtracted  from  the  current 
measurement,  resulting  in  the  state  estimation  error  commonly  called  the  innovation  or  innovation 
process.  The  innovation  is  then  “weighted”  by  the  Kalman  gain  matrix  (K),  which  is  calculated  in  the 
Kalman  filter  covariance  loop  shown  in  Figure  3-6. 


Figure  3-6  Kalman  Filter  Covariance  and  Kalman  Gain  Loop 


As  stated  above,  the  Kalman  gain  is  calculated  in  the  covariance  loop  of  the  Kalman  filter  shown 
in  Figure  3-6.  It  should  be  emphasized  that  the  Kalman  gain  is  independent  of  the  state.  This  can  be  seen 
by  noticing  that  the  only  inputs  to  the  covariance  loop  are  the  process  (or  plant)  error  covariance  (Q),  the 
measurement  (or  observation)  error  covariance  (R),  and  its  previously  calculated  Kalman  gain. 
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Figure  3-7  Chi-Square  Gate  Test  Block  Diagram 


Figure  3-7  illustrates  the  internal  block  diagram  of  the  chi-squared  gate  test  block  from  Figure  3- 
4.  This  block  performs  the  correlation  criterion  calculations  described  in  equation  (2.24)  and  compares 
this  result  to  the  gate  size.  Looking  forward  to  the  analysis  of  results  section  of  this  paper,  the  plots 
concerning  the  gating  criterion  correspond  to  d2  in  Figure  3-7.  The  MatLab®  function  that  performs  the 
inverse  operation  here  is  the  pseudo-inverse  function  with  an  imposed  tolerance  of  10'6.  The  pseudo¬ 
inverse  function  is  used  instead  of  the  inverse  or  division  functions  due  to  the  possibility  of  singularity  of 
the  innovation  covariance  matrix  which  may  result  from  round-off  errors. 

The  matrix  multiplication  blocks  in  Figure  3-7  require  knowledge  of  the  size  of  the  matrices  to  be 
multiplied.  These  values  cannot  be  changed  during  a  simulation;  this  is  the  source  of  difficulty  in  passing 
varying  sized  matrices  that  would  be  required  for  representing  missed  measurements  or  extraneous  clutter 
(false  alarm)  measurements.  One  may  be  able  to  define  an  initial  arbitrarily  large  matrix  capable  of 
containing  the  maximum  probable  number  of  targets  measured  at  any  instant.  In  the  event  of  a  missed 


30 


target,  the  measurement  value  could  be  represented  by  not  a  number  (NaN).  The  algorithms  that  operate 
on  the  measurements  would  then  need  to  be  capable  of  detecting  the  NaN  and  either  suspend  calculation, 
or  hold  the  previous  calculated  value. 

3.3  Assumptions 

A  number  of  assumptions  are  asserted  in  order  to  simplify  the  creation  of  the  simulation  model. 
These  assumptions  may  be  grouped  into  the  following  topics: 

3.3.1  Fixed  Sample  (or  Measurement)  Intervals 

Although  a  fixed  measurement  interval  is  assumed,  an  actual  rotating  radar  antenna  generally  will 
not  measure  a  moving  target  at  a  fixed  interval  unless  the  target  path  is  strictly  radial.  This  assumption  is 
made  primarily  due  to  the  difficulty  of  implementing  varying  sample  times  in  Simulink®.  It  should  be 
noted  however,  the  Kalman  filter  algorithm  doesn’t  required  fixed  sample  intervals,  thus  future 
amendments  to  this  model  may  remove  this  assumption. 

3.3.2  Single  Target 

A  single  target  is  assumed  in  this  analysis.  This  is  done  since,  as  eluded  to  earlier,  each  target 
requires  a  dedicated  Kalman  filter.  Additionally,  the  influences  of  biases  in  the  single  target  case  will 
also  apply  for  multiple  targets  assuming  that  dedicated  Kalman  filters  are  required  for  each  target.  This 
assumption  does  however,  exclude  the  analysis  of  situations  such  as  crossing  targets,  and  closely 
separated  targets  (formations). 

3.3.3  No  Clutter  /  False  Alarms 

The  assumption  that  no  false  alarms  occur  is  made  for  reasons  similar  to  the  single  target 
argument.  In  the  event  that  clutter  generates  false  alarms,  a  bank  of  Kalman  filters  is  needed  to  process 
all  potential  tracks  since  it  cannot  be  determined  if  the  detection  is  due  to  clutter  or  a  real  target.  Future 
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amendments  to  this  model  would  gain  from  the  vectorization  of  the  Kalman  filter  algorithm.  This  may  be 
accomplished  by  implementing  an  S-function  capable  of  processing  varying  width  vectors. 

3.3.4  100  %  Probability  of  Detection 

It  is  assumed  that  each  target  is  detected  once  on  each  scan.  While  this  assumption  rarely  reflects 
reality,  the  difficulty  of  processing  missed  detections  necessitated  its  acceptance.  A  missed  detection 
must  be  represented  mathematically  in  the  simulation  model.  At  each  observation,  the  Kalman  filter 
accepts  the  target’s  spatial  coordinates  and  performs  the  ensuing  calculations.  Representing  the  missed 
detection  as  a  vector  of  Os  impacts  the  calculation  of  the  predicted  state  to  varying  degrees,  depending  on 
the  previous  measurement  and  state.  To  remove  this  assumption,  some  method  must  be  used  to  alleviate 
this  problem. 
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4.  Analysis  of  Results 


4.1  Overview 


x  ar  «  1 

Figure  4-1  Representative  Constant-Error  Ellipses 


Refer  to  Figures  4-1  and  4-2  to  help  visualize  the  following  discussion.  A  substantia]  portion  of 
current  literature  assumes  that  the  cross-range  semi-axis  (ra9)  is  much  greater  than  the  down-range  semi¬ 
axis  (cxr)  when  analyzing  worst  case  radar-target  geometries[Bar90].  In  an  effort  to  model  a  realistic 
networked  surveillance  radar  system,  the  following  measurement  statistics  are  chosen:  a9  =  0.18  degrees 
(deg)  =  Tt-10'3  radians  (rads)  and  ar  =  0.125  nautical  miles  (nmi)  =  231.5  meters  (m)[Bar90],  Recalling 
that  the  semi-axes  of  the  constant-error  ellipses  are  equal  to  ar  and  r-cr9,  and  considering  the  self-imposed 
maximum  detectable  range  of  approximately  100  Km,  this  corresponds  to  r-c9  =  100,000  m  •  7I-10'3  rads  = 
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314  m.  Due  to  the  target  path  geometries  chosen,  the  minimum  down-range  distance  is  nearly  35,000  m, 
resulting  in  r  ae  =  35,000  m  •  tc-  1  O'3  rads  =  109  m. 


Figure  4-2  Constant-Error  Ellipses  vs.  Target  Position 


The  down-range  measurement  variance  is  not  a  function  of  range  but  of  the  specific  parameters  of 
the  radar  such  as  pulse  width,  number  of  pulses  integrated,  Doppler  frequency,  clock  stability,  etc. 

Because  of  this  fact,  the  range  independent  ar  =231 .5  m.  The  range  at  which  the  measurement  error 
ellipse  becomes  circular  is  r  =  ar/a9  =  231.5/(71-1  O'3)  =  73,688  m. 


Having  determined  the  constant-error  ellipse  dimensions,  it  is  observed  that  the  ellipses  will  not 
be  highly  eccentric,  but  rather  oval  or  circular  contours.  Furthermore,  the  area  (A)  of  the  constant-error 
ellipses  can  be  calculated  by  A  =  7t-r-ae-ar.  It  is  obvious  that  the  maximum  area  (Amax)  will  occur  at  the 
maximum  detectable  range,  Amax  =  n  •  100,000  m  •  re- 1 0‘3  rads  •  231.5  m  =  228,481  m2.  Likewise,  the 
minimum  area  (Amin)  will  occur  at  a  down-range  distance  of  35,000  m,  then  Amm  =  7t  ■  35,000  m  •  Tt-10'3 
rads  •  231.5  m  s  80,000  m2. 

All  of  the  following  simulations  are  performed  thirty  times  for  each  increment  of  the 
independent  variable,  the  minimum  required  to  assert  the  Gaussian  assumption  in  the  confidence  level 
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calculations.  The  initial  seed  for  each  random  number  generator  is  uniquely  determined  for  each  run  as  a 
function  of  the  system  clock. 


4.2  Positional  Bias 

Radars  A  and  B  are  positioned  at  (50  Km,  100  Km)  and  (50  Km,  0  Km)  in  global  Cartesian 
coordinates,  respectively,  and  the  global  origin  is  at  (0,  0),  as  shown  in  Figure  4-3.  The  target  travels  in  a 
diagonal  path  at  200  m/sec  in  the  -x  and  -y  directions.  The  initial  detection  of  the  target  is  at  (100  Km, 
100  Km)  in  global  Cartesian  coordinates.  Thus  the  target  is  centered  between  radars  A  and  B  after  250 


(100  Km,  100  Km) 


(50  Km,  100  Km)  t  =  0sec. 

Y  Vel.  =  (-200m/s,  -200m/s) 


(-20  Km,  -20  Km) 
t  =  600  sec 


Figure  4-3  Position  Bias  Target  Path 


seconds  have  elapsed  after  initial  detection.  The  injected  biases  at  either  of  the  radars  A  or  B  are  such 
that  they  are  directed  along  the  x-axis.  The  measurement  error  ellipses  at  this  point  are  greater  in  the 
down-range  direction  than  they  are  in  the  cross-range  direction.  In  the  absence  of  the  cross-range  biases, 
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these  ellipses  “line-up”,  however  due  to  the  biases,  they  are  misaligned  by  relative  value  of  the  bias 
between  radars  A  and  B. 

When  calculating  the  upper  bound  on  the  position  bias,  the  minimum  detected  range  is  about  35 
Km,  however  the  range  occurring  during  the  worst-case  geometry  is  50  Km.  If  the  smaller  down-range 
distance  is  used,  jbp|=l  .34-35  Km-7t-10'3  =  147  m.  Alternatively  using  the  worst-case  range  of  50  Km, 
|bp|=l  .34-50  Km-7i-10‘3  =  210  m.  As  illustrated  in  Figure  4-4,  the  simulated  |bp|  =  180  m  falls 
approximately  at  the  midpoint  of  these  two  values. 

Since  these  positional  biases  cause  the  same  misalignment  everywhere  in  the  detectable  area,  the 
worst-case  geometry  is  determined  merely  by  the  eccentricity  of  the  ellipses  and  the  orientation  of  the 
ellipses  with  respect  to  the  biases. 

It  can  be  seen  in  Figure  4-4  that,  as  predicted  by  our  previous  development  of  Chapter  two,  the 
worst-case  radar-target  geometry  occurs  when  the  target  is  centered  between  the  two  radars  and  the 
positional  bias  is  perpendicular  to  the  line  joining  the  two  radars.  It  should  be  noted  that  position  biases 
do  not  generate  as  pronounced  effects  in  the  worst-case  geometry  as  do  the  range  and  azimuth  biases,  as 
will  be  shown.  Lambda,  the  variable  along  the  vertical  axis  in  Figure  4-4,  represents  the  total  normalized 
bias,  or  equivalently,  the  noncentrality  parameter  in  the  Chi-square  distribution. 

The  transient  properties  of  the  Kalman  filter  cause  the  small  initial  Chi-square  test  results.  This 
occurs  because  the  initial  confidence  in  the  measurements  is  small,  so  the  measurements  (and  added 
position  bias)  are  weighed  less  heavily.  In  this  case,  the  filter  relies  more  upon  its  predicted  measurement 
which  is  derived  from  the  filter’s  linear  target  path  model.  It  should  be  noted  here  that  the  simulation 
results  depend  to  some  extent  on  the  initial  values  of  the  prediction  error  covariance  matrix,  P0,  and  the 
initial  state  values,  S0.  The  a  priori  knowledge  of  the  initial  target  location  and  velocity  allows  determines 
the  S0,  and  the  process  and  measurement  statistics  determine  P0. 
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Figure  4-5  Position  Bias  95%  Confidence  Level  (one-sided)  vs.  Target  Position 
4.3  Range  Bias 

Figure  4-6  illustrates  the  radar-target  geometry  chosen  for  the  range  bias  simulation.  This 
geometry  is  chosen  to  minimize  the  initial  convergence  effects  of  the  filter,  and  to  allow  the  target  path  to 
pass  through  the  (50  Km,  50  Km)  position  to  observe  the  effects  of  an  orthogonal  radar-target  geometry. 

As  with  the  positional  biases,  the  worst-case  geometry  for  range  biases  occurs  for  a  target 
centered  between  the  two  radars.  The  biases  in  this  case  however,  are  down  range  biases  of  equal  value 
resulting  in  a  misalignment  of  the  measurement  error  ellipses  in  varying  directions  depending  on  target 
position.  The  theoretical  range  bias  required  at  each  radar  to  yield  the  non-centrality  parameter  X  =  1 .2  is, 
|br|=0.67  •  CTr  =  1 55  m.  Figure  4-7  illustrates  the  simulated  bias  value  necessary  to  reach  the  non-centrality 
value  is  somewhat  higher,  nearly  1 80  m.  This  may  be  explained  by  the  fact  that  the  range  bias  shifts  the 
error  ellipses  along  their  major  semi-axes  since  the  down-range  variance  is  greater  than  the  cross-range 
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error  in  this  position.  Meanwhile,  the  theoretical  bound  on  the  range  bias  assumes  the  cross-range  error  is 
much  greater  than  the  down-range  variance.  At  the  worst-case  position,  ar  =  23 1 .5  m,  while  r-ae  =  35  Km 
■ti-10'3  =  110  m. 


(100  Km,  100  Km) 
t  =  0  sec. 

Y  Vel.  =  (-200m/s,  -200m/s) 


(“20  Km,  “20  Km) 
t  =  600  sec 


Figure  4-6  Range  Bias  Target  Path 
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4.4  Azimuth  Bias 


For  the  analysis  of  the  azimuth  bias  effects,  the  radar-target  geometry  illustrated  in  Figure  4-3  is 
again  chosen  to  place  the  worst-case  target  position  near  the  center  of  the  path.  As  discussed  in  Chapter 
two,  since  the  down-range  variance  is  greater  than  the  cross-range  variance,  the  worst-case  case  radar- 
target  geometry  occurs  at  the  midpoint  of  the  line  joining  the  two  radars.  The  effect  of  the  azimuth  biases 
at  this  worst-case  position  is  similar  to  the  effect  of  position  bias,  since  the  measurement  will  be  shifted,  if 
not  along  the  x-axis,  at  least  tangential  to  the  x-axis.  The  bound  for  the  allowable  azimuth  bias  is  then 
calculated  as  |be|=0.67-ae  =  0.67-0.18  deg  =  0.12  deg. 

Figure  4-9  reveals  the  azimuth  bound  to  be  approximately  0.1 1  deg,  which  seems  reasonable  in 
light  of  Figure  4-10  which  shows  that  the  one-sided  95%  confidence  interval  is  nearly  0.17.. 


Figure  4-9  Azimuth  Bias  vs.  Target  Position 
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Figure  4-10  95%  Confidence  Level  of  Azimuth  Bias  vs.  Target  Position 

4.5  Simultaneous  Position-Range- Azimuth  Biases 

Finally,  the  radar-target  geometry  shown  in  Figure  4-6  is  chosen  for  the  simultaneous  position, 
range,  and  azimuth  bias  simulation  such  that  the  influence  of  an  orthogonal  radar-target  geometry  may  be 
examined,  if  they  exist.  The  worst-case  radar-target  geometry  for  all  three  types  of  biases  happens  to 
occur  at  the  same  position,  namely  at  the  midpoint  of  the  line  joining  the  two  radars.  Recall  that  the 
position  bias  is  a  relative  bias,  effectively  shifting  the  measurements  an  equal  amount  regardless  the 
position  of  the  target.  Contrary  to  the  position  bias,  the  range  and  azimuth  biases  affect  the  measurements 
differently  depending  on  the  target’s  position  relative  to  the  radar. 

Since  the  azimuth  bias  seems  to  have  the  most  position-dependent  properties,  the  simultaneous 
effect  of  the  three  biases  is  analyzed  by  setting  the  position  and  range  biases  at  their  upper  bound,  then 
assign  the  azimuth  bias  as  the  independent  variable. 
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That  is  |bp|  =  0.77-r-Oe  =  0.77  •  35  Km-7t-10'3  =  85  m,  since  the  worst-case  position  coincides  with 
the  minimum  detectable  range  for  this  flightpath.  The  range  bias  is  calculated  as  |br|  =  0.39-or  = 

0.39-23 1 .5  m  =  90  m.  The  upper  limit  on  the  azimuth  bias  is  calculated  as  |be|  =  0.39-ae  =  0.39-0. 1 8  deg  = 
0.07  deg. 

The  results  of  the  simultaneous  position,  range,  and  azimuth  simulation  illustrated  in  Figure  4-1 1, 
compares  favorably  with  the  theoretical  simultaneous  upper  bounds  on  each  of  the  biases  in  consideration 
of  the  one-sided  confidence  interval  of  Figure  4-12.  The  simulated  azimuth  bias  required  to  exceed  the 
upper  bound  on  the  correlation  criterion  is  |be|  =  0.06  deg.  It  may  be  observed  that  the  result  of 
combining  these  biases  tends  to  “blend”  the  spatial  properties  of  each,  that  is  the  relatively  position- 
independent  property  of  the  position  bias  and  the  more  position-dependent  properties  of  the  range  and 
azimuth  combine  to  form  a  moderately  position-dependent  response.  Of  course  the  convergence  of  the 
Kalman  filter’s  covariance  matrix  causes  the  “slow”  transient  response  immediately  after  initial  detection. 
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Figure  4-11  Simultaneous  Position-Range- Azimuth  Bias  vs.  Target  Position 


Time  (sec)  Azimuth  Bias  (deg) 

Figure  4-12  Position-Range-Azimuth  95%  Confidence  Level 
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5.  Summary  /  Recommendations 


This  thesis  examines  the  effect  of  registration  errors,  namely,  position,  range,  and  azimuth  biases 
both  individually  and  collectively.  The  worst-case  radar-target  geometries  may  be  predicted  given 
knowledge  of  the  down-range  and  cross-range  measurement  error  statistics.  The  quantitative  effects  of 
the  biases  may  be  analyzed  by  employing  statistical  correlation  gates.  These  elliptical  correlation  gates 
are  defined  by  the  Chi-square  statistics  of  the  measurement  errors.  The  degree  of  ellipticity  of  these 
equal-probability  contours  are  governed  by  the  measurement  error  statistics  and  the  down-range  distance. 

5.1  Coordinate  Transformations 

To  gain  an  advantage  from  measurements  at  multiple  distributed  radars,  the  measurements  must 
be  transformed  into  a  common  reference  frame.  To  remove  the  inherent  weakness  of  two-dimensional 
radar  to  determine  correct  slant  range,  target  returns  are  generated  in  three  dimensions.  The  individual 
radar  measurements  are  then  transformed  into  “global”  coordinates.  First  the  radar  geographical 
coordinates  are  transformed  to  geocentric  coordinates  using  an  ellipsoid  to  approximate  true  earth  geoid 
shape.  The  local  coordinate  system  of  one  radar  is  then  transformed  to  that  of  another,  or  alternatively  to 
an  arbitrary  “global”  coordinate  system  through  a  translation  and  a  series  of  rotations,  (see  Appendix  3). 

5.2  Kalman  Algorithm 

Presented  in  Appendices  1  and  2  is  the  development  of  the  discrete  Kalman  filter  difference 
equations  which  provide  the  minimum  variance,  or  most  accurate  estimate  of  the  actual  target  position. 
The  Kalman  filter  employs  an  embedded  model  of  the  system,  which  in  this  case  is  a  target  traversing  a 
linear  path  in  cartesian  coordinates.  The  Kalman  filter  is  a  recursive  estimator,  with  the  advantage  that 
the  required  memory,  or  sample  length,  does  not  increase  with  time.  All  the  needed  information  is  stored 
in  the  state  and  error  covariance  matrices. 
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Two  approaches  to  the  development  of  the  Kalman  filter  equations  are  provided.  The  first 
assumes  that  the  noise  processes  and  the  initial  state  are  normally  distributed,  and  the  mean-square 
estimate  criterion  is  then  used  to  obtain  an  optimum  filter  which  estimates  the  state  conditioned  on  the 
measurement.  From  the  assumptions  of  a  linear  model  and  Gaussian  statistics,  a  linear  filter  results. 

The  second  approach  is  the  linear  mean-square  estimator,  which  makes  no  assumptions  about  the 
process  distribution  functions.  This  approach  makes  use  of  the  orthogonality  between  the  estimation  error 
and  the  observations,  the  concept  of  the  innovation  sequence,  and  a  recursive  updating  formula  for  the 
estimate  when  a  new  measurement  is  available. 

5.3  Analytical  Model  Development 

Given  that  the  radar  measurements  are  independent  Gaussian  distributed  random  variables,  the 
mean-square  error  criterion  is  a  Chi-square  distributed  random  variable,  i.e.,  the  sum  of  squared 
independent  Gaussian  random  variables  is  a  Chi-square  distributed  random  variable.  The  maximum 
allowable  bias  is  determined  by  choosing  a  minimum  probability  of  correlation  between  two  radar 
measurements  based  on  the  desired  system  track  accuracy.  Each  imposed  bias  is  then  viewed  as  the  non¬ 
centrality  parameter  of  a  non-central  Chi-square  distributed  random  variable,  which  allows  for  the 
determination  of  whether  or  not  the  two  measurements  “probably”  come  from  the  same  target.  Given  this 
metric,  the  allowed  biases  in  position,  range,  and  azimuth  may  be  quantitatively  determined. 

5.4  System  Model 

Given  that  tracking  is  performed  for  a  target  with  a  linear  path  (underlying  assumption),  the 
appropriate  state  transition  matrix  is  a  linear  function  of  the  previous  position  and  the  product  of  target 
velocity  and  observation  interval.  The  radars  measure  only  position,  however,  to  make  a  more  accurate 
prediction  the  state  vector  is  augmented  to  include  velocities  in  the  x  and  y  directions.  The  observation 
matrix,  therefore,  merely  extracts  the  position  components  of  the  state  vector. 
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Due  to  the  linearity  of  the  system,  the  next  state  is  determined  by  the  superposition  of  the  effect 
of  the  actual  state  value  and  the  effect  of  the  input  noise.  The  state-transition  matrix  accounts  for  the 
contribution  to  the  predicted  state  of  the  past  history  of  the  system  “stored”  in  the  present  state.  The 
observations  are  linear  combinations  of  the  state  components,  which  are  corrupted  by  additive  noise. 

5.5  Bias  Effects 

As  predicted  by  the  analytical  model,  position  bias  is  relative.  That  is,  if  both  radars  have  the 
same  vector  bias,  their  measurements  fall  within  the  correlation  gate  most  of  the  time,  no  matter  how 
large  the  bias.  Thus  for  the  worst-case  bias,  antipodal  biases  are  imposed,  i.e.,  opposing  vector  biases,  for 
each  of  the  two  radars.  Notice  that  the  position  biases  impose  a  relatively  constant  correlation  criterion 
regardless  of  the  relative  radar-target  geometries.  The  nature  of  the  constant-error  contours  determines 
the  worst-case  radar-target  geometries  and  the  vector  direction  of  the  worst-case  biases.  The  range 
variance  is  larger  than  the  cross-range  variance  for  the  chosen  measurement  statistics  and  radar  coverage 
area.  Thus,  worst-case  positional  bias  occurs  when  the  target  is  centered  between  radars  and  the  bias  is 
perpendicular  to  the  imaginary  line  joining  the  two  radars. 

The  worst-case  geometry  for  range  biases  is  the  same  as  for  position  biases,  however,  range 
biases  are  not  relative,  and  each  radar  can  contribute  an  error.  Whether  these  errors  are  additive  or 
partially  cancel  each  other  depends  on  the  exact  geometry  and  on  the  sign  of  the  bias,  which  explains  why 
range  bias  shows  more  position  dependence  than  does  position  bias. 

The  worst-case  geometry  for  the  azimuth  biases  coincides  with  the  position  and  range  worst-case 
geometries.  The  azimuth  biases  are  also  more  position  dependent  than  the  position  biases.  In  the 
imposed  radar-target  geometry,  the  cross-range  variance  is  minimum  when  the  target  is  at  the  mid-point 
of  the  line  connecting  the  two  radars.  Thus,  the  worst-case  bias  scenario  is  for  both  radars  to  have 
azimuth  biases  in  the  same  rotational  direction. 
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Combining  all  three  biases  simultaneously  results  in  a  more  position  dependent  response.  Of 
course  in  the  general  case  the  effect  of  simultaneous  position,  range,  and  azimuth  biases  varies  depending 
upon  radar  measurement  statistics,  down-range  distance,  radar-target  geometries,  and  the  relative  vector 
directions  of  the  biases. 

5.6  Recommendations 

The  implemented  model  allows  for  the  simulation  of  one  target  tracked  by  two  spatially 
distributed  radars,  yet  much  may  be  done  to  improve  the  usefulness  of  this  model.  The  following 
paragraphs  present  a  number  of  recommendations  and  discusses  their  possible  implementation. 

5. 6. 1  Vectorize  the  Kalman  filter  algorithm 

To  process  all  incoming  signals  from  the  radar  detection  hardware,  not  one  Kalman  filter  but  an 
entire  bank  of  filters  must  be  used  to  formulate  the  statistics  for  all  potential  targets,  whether  they  are 
from  actual  targets  or  from  false  clutter  detections.  Some  functions  in  MatLab®  accept  variable  vector 
“widths”  as  input.  It  would  be  advantageous  to  “vectorize”  a  Kalman  algorithm  function  so  that  at  each 
discrete  moment  any  detection  within  the  radar  range-azimuth  detection  gate  could  be  processed  as  a 
confirmed  track,  a  tentative  track,  or  a  new  detection. 

The  vectorization  may  be  realized  by  representing  the  Kalman  filter  difference  equations  in  their 
z-transform  transfer  function  form.  The  inclusion  of  the  initial  state  values  would  be  necessary  for  this 
method  to  correctly  represent  the  Kalman  filter  algorithm.  Additionally,  the  measurement  statistics  must 
be  observable  at  each  instant  to  allow  execution  of  the  Chi-square  gate  test. 

5. 6.2  Accept  ASCII/MA  T  plot  files 

As  stated  earlier,  output  plot  (or  detection)  files  from  accepted  radar  simulation  software 
packages  should  be  acceptable  input  to  the  simulation  model.  In  MatLab®,  .mat  files  are  ASCII  text  files 
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formatted  such  that  the  monotonic  increasing  or  decreasing  independent  variable  is  in  the  first  row  of  the 
file  and  the  accompanying  dependent  variables  are  in  the  remaining  rows.  In  the  event  that  the  input  is 
already  in  this  format,  integration  is  easily  executed  with  the  load  command.  However,  if  the  data 
provided  is  in  another  format,  the  primary  difficulty  is  converting  the  file  to  .mat  format. 

5.6.3  Multiple  Radars 

In  the  event  that  more  than  two  radars  detect  a  target  within  the  same  correlation  gate  in  a 
networked  radar  system,  the  measurement  statistics  from  all  radars  may  be  processed  together  to  take  full 
advantage  of  sensor  fusion.  This  joint  processing  could  be  performed  by  a  series  of  cascaded  Kalman 
filters,  in  which  the  first  stage  is  a  bank  of  Kalman  filters  and  each  filter  processes,  (fuses),  all  possible 
combinations  of  detected  targets  in  a  correlation  gate.  The  number  of  filters  required  in  the  ensuing  levels 
of  the  cascade  would  of  course  decrease  by  a  factor  of  two.  Alternatively,  after  the  first  level  of  the 
cascade  the  statistical  distances  may  be  compared  to  select  the  smallest  distance  so  that  the  track  may  be 
updated  with  a  filtered  estimate. 

5. 6.4  Degradation  of  Tracking  Accuracy 

This  thesis  addresses  the  worst  case  situation  in  a  networked  radar  system,  i.e.,  the  case  where 
biases  force  the  initiation  of  additional  tracks  for  a  single  target.  As  discussed  in  the  introduction,  biases 
may  introduce  a  continuum  of  tracking  degradation  effects,  from  reduced  accuracy  relative  to  maximum 
theoretical  accuracy  to  the  creation  of  a  separate  track.  To  actually  estimate  the  biases  of  each  radar  and 
correct  for  them  is  a  much  more  difficult  problem.  It  may  be  shown  that  the  ability  to  accurately  estimate 
these  biases  is  proportional  to  the  number  of  available  targets  in  overlapping  coverage  areas  and  the  radar- 
target  geometries  [Bar90], 

The  Munkres  algorithm  mentioned  earlier,  provides  the  optimum  assignment  in  the  event  of  N 
tracks  and  M  detections,  (where  there  may  be  more  or  less  tracks  than  detections).  The  Munkres 
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algorithm  is  essentially  a  recursive  matrix  reduction  algorithm  which  performs  subtraction  in  row  or 
column  space  while  continually  testing  for  the  “closest”  detection  for  each  track[Bla86].  Although  the 
results  reported  here  do  not  use  the  Munkres  algorithm,  MatLab®  code  implementing  this  algorithm  is 
included  in  Appendix  4. 

5.6.5  The  Effect  of  Sample  Time 

It  may  be  informative  to  investigate  the  behavior  of  error  covariances  when  the  observation 
interval  becomes  longer.  More  precisely,  it  may  be  worth  determining  whether  finite  limit  values  exist 
for  the  error  covariances  as  the  observation  interval  becomes  large  and  under  what  conditions  they  are 
independent  of  the  initial  value  of  the  error  covariance. 

5.6.6  Multiple-Hypothesis  Tracking 

Although  a  constant  velocity  target  is  assumed  here,  many  situations  may  occur  where  the  target 
maneuvers  significantly  for  a  period  of  time.  In  this  case  the  assumed  linear  system  model  is  no  longer 
valid,  instead  a  constant  acceleration  model  may  be  implemented  when  the  target  no  longer  falls  within 
the  correlation  gate.  When,  (or  if),  the  target  resumes  a  constant  velocity  profile,  the  assumed  model  may 
revert  back  to  the  linear  model.  In  essence  then,  it  can  be  assumed  or  hypothesized  that  the  target  path 
has  multiple  “modes”  or  profiles,  hence  the  name  multiple-hypothesis  tracking.  Of  course  multiple 
modes  are  more  difficult  to  implement  than  the  single  system  model,  and  problems  linked  to  this 
switching  process  must  be  solved  to  maintain  tracking  continuity  between  multiple  system  models. 
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1A.  Appendix  1  Kalman  Filtering  Introduction 


1A.  1  Linear  System  Model 

In  many  estimation  problems,  an  unknown  vector  represents  the  time  evolution  of  a  system,  and 
the  measurements  Z  are  the  observed  output  from  the  system.  The  state-variable  approach  is  a  valuable 
method  to  describe  such  dynamic  systems  where  the  input-output  relationship  is  described  in  the  time 
domain  by  a  state-transition  model  together  with  an  output  observation  model.  The  state  represents  the 
internal  condition  of  the  system  and  accounts  for  its  memory  of  past  inputs.  The  inputs  may  consist  of 
deterministic  functions  of  time  together  with  stochastic  processes  representing  unpredictable  variables  or 
noise.  The  output  is  a  function  of  the  state,  usually  corrupted  by  random  measurement  errors. 

A  discrete-time  system  is  described  by  the  following  pair  of  equations: 

Sk+1=F(Sk,Uk,  Vk,  k)  (1A.1) 

Zk+i  =  G(Sk+1,  Wk+i,  k+1)  (1A.2) 

Where  the  integer  subscripts  (k,k+l)  denote  discrete-time  instants, 

Sk  is  an  n-dimensional  vector  representing  the  state  at  time  k, 

Uk  is  a  p-dimensional  vector  representing  the  deterministic  control  inputs, 

Vk  is  the  q-dimensional  random  vector  representing  the  model  process  (or  plant)  noise, 

Wk  is  an  r-dimensional  random  vector  representing  the  observation  (or  measurement)  noise,  and 
Zk  is  the  m-dimensional  vector  of  the  observations  (or  measurements). 
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F(.)  accounts  for  the  state  transition  from  Sk  to  Sk+1  due  to  the  input  (Uk,  Vk)  occurring  at  time  k.  G(.) 
performs  the  projection  of  the  internal  state  on  the  output  measurable  variables  Zk+1,  as  well  as  the  effect 
of  the  random  errors  Wk+].  The  explicit  dependence  on  k  accounts  for  a  possibly  non-stationary  system. 

A  fundamental  case  to  be  considered  is  that  of  linear  systems,  primarily  because  of  their  practical 
interest  and  also  because  complete  results  in  estimation  theory  are  available  for  this  class  of  systems. 
Equations  (1A.1)  and  (1A.2)  particularize  to  the  following: 

Sk+i  =  ®kSk  +  BkUk  +  GkVk  (1A.3) 

Zk+,  =  Hk+]Sk+i  +  Lk+1Wk+1  (1A.4) 

where  Ok,  Bk,  Gk,  Hk+i,  Lk+)  are  real-valued  matrices  having  dimensions  (nxn),  (nxp),  (nxq),  (mxn), 

(mxr),  respectively.  Equation  (1A.3)  shows  that,  because  of  the  linearity  of  the  system,  the  next  state  Sk+i 
is  determined  by  the  superposition  of  two  terms:  the  effect  of  the  actual  state  value  Sk  and  that  of  the  input 
samples  Uk  and  Vk.  The  state-transition  matrix  Ok,  accounts  for  the  contribution  to  Sk+i  of  the  past  history 
of  the  system,  stored  in  Sk.  In  particular,  when  input-free  conditions  occur  (i.e.  Uk  =  0  and  Vk  =  0),  the 
state  depends  on  the  initial  condition  So  only,  according  to  the  equation 

Sk+1  =  3>k  •  <I>k.i . 0|  •  <t>o  •  So  (1A.5) 

In  this  case,  the  stability  of  the  system,  i.e.  the  condition  for  Sk  to  remain  bounded,  depends  on  the 
matrices  O,  (j  =  0,1, — ,k).  The  matrices  Bk  and  Gk  account  for  the  possibility  of  the  system  being  forced 
into  a  new  state  by  input  variables.  Equation  (1A.4)  shows  that  the  observations  Zk+I  are  linear 
combinations  of  the  state  components,  corrupted  by  additive  noise.  The  matrix  Bk+1  indicates  how  the 
components  of  Sk+]  are  combined  to  form  the  observed  vector  Sk+]  and  therefore  accounts  for  the  ability 
of  the  system  state  to  be  determined  from  the  measurements. 
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In  order  to  complete  the  description  of  discrete-time  linear  systems,  it  is  worth  deriving  the  time 
evolution  of  the  statistics  of  Sk  and  Zk  when  a  random  input  Vk  is  present.  Let  Vk  and  Wk  be  zero-mean, 
white  processes  having  the  following  covariance  matrices: 


E{VjVkT}  =  Qk§k, 

(1A.6) 

E{WjWkT}  =  Rk5kj 

(1A.7) 

E{VjWkT}=0 

(1A.8) 

where  5kj  is  the  Kroenecker  operator  (5kj  =  0  if  k  *  j  and  5^=1)  while  Qk  and  Rk  are  positive  semi- 
definite  real  matrices.  These  matrices  describe  the  correlation  between  the  different  components  of  Vk 
and  Wk,  at  the  same  time  sample  k,  whereas  no  correlation  exists  between  samples  taken  at  different 
times  (white  processes).  The  expected  value  of  Sk  is  simply  obtained  from  equation  (1  A. 3)  as 

E{Sk+1}=Ok-E{Sk}+Bk-Uk  (1A.9) 

And  similarly  for  Zk, 

E{Zk}=Hk-E{Sk}  (1A.10) 

The  covariance  matrix  of  Sk  (representing  the  correlation  between  the  components  of  Sk  at  the  same 
instant  k)  is  obtained  from  equation  (1  A.3)  and  satisfies  the  following  recursive  relation: 

Cov{Sk+i}  =  Ok-cov{Sk}-OkT  +  Gk-Qk-GkT  (1A.11) 

Similarly,  the  covariance  of  Zk  is  given  by 

Cov{Zk}  =  Hk-cov{Sk}-HkT  +  Lk-Rk-LkT  (1  A.  12) 

1A.2  Discrete-Time  Kalman  Filtering 
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Before  describing  the  Kalman  filtering  algorithm,  it  may  help  to  consider  an  intuitive  example  of 
recursive  estimation[Bro92],  Consider  the  problem  of  estimating  the  mean  of  some  random  constant 
based  on  a  sequence  of  noisy  measurements.  Let  us  assume  that  our  estimate  is  to  be  the  sample  mean 
and  that  we  wish  to  refine  our  estimate  with  each  new  measurement  as  it  becomes  available.  That  is, 
think  of  processing  the  data  on-line.  Let  the  measurement  sequence  be  denoted  as  z\,  z2, . . . ,  z„,  where 
the  subscript  denotes  the  time  at  which  the  measurement  is  taken.  One  method  of  processing  the  data 
would  be  to  store  each  measurement  as  it  becomes  available  and  then  compute  the  sample  mean  in 
accordance  with  the  following  algorithm: 

First  measurement  zt:  Store  z\  and  estimate  the  mean  as 


ml=z] 


(1A.13) 


Second  measurement  z2:  Store  z2  along  with  Z\  and  estimate  the  mean  as 


Zj  “f“  Z2 


m2  j 


(1A.14) 


Third  measurement  z3:  Store  z3  along  with  Z\  and  z2  and  estimate  the  mean  as 


m,  = 


Zi+Zo  +  Zt 


3  -  3 


(1A.15) 


And  so  forth. 

Clearly,  this  would  yield  the  correct  sequence  of  sample  means  as  the  experiment  progresses.  It  should 
also  be  clear  that  the  amount  of  memory  needed  to  store  the  measurements  keeps  increasing  with  time, 
and  also  the  number  of  arithmetic  operations  needed  to  form  the  estimate  increases  correspondingly.  This 
would  lead  to  obvious  problems  when  the  total  amount  of  data  is  large.  Thus,  consider  a  simple  variation 
in  the  computational  procedure  in  which  each  new  estimate  is  formed  as  a  blend  of  the  old  estimate  and 
the  current  measurement.  To  be  specific,  consider  the  following  algorithm; 
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First  measurement  z,:  Compute  the  estimate  as 


mx  =  z, 


(1A.16) 


Store  m]  and  discard  z. . 

Second  measurement  z2:  Compute  the  estimate  as  a  weighted  sum  of  the  previous  estimate  mx  and  the 
current  measurement  z2: 


m2  =\mx  +jz2 


(1A.17) 


Store  m2  and  discard  z2  and  m] . 

Third  measurement  z3:  Compute  estimate  as  a  weighted  sum  of  m2  and  z3: 

W3=jOT2+iZ3  (1A.18) 


Store  m3  and  discard  z3  and  m2 . 

And  so  forth.  It  should  be  obvious  that  at  the  k,h  state  the  weighted  sum  is 

% =(¥j"vi+(i-K  <1A-19) 

Clearly,  the  above  procedure  yields  the  same  identical  sequence  of  estimates  as  before,  but  without  the 
need  to  store  all  the  previous  measurements.  We  simply  use  the  result  of  the  previous  step  to  help  obtain 
the  estimate  at  the  current  step  of  the  process.  In  this  way,  the  previous  computational  effort  is  used  to 
good  advantage  and  not  wasted.  The  second  algorithm  can  proceed  on  ad  infinitum  without  a  growing 
memory  problem.  Eventually,  of  course,  as  k  becomes  extremely  large,  a  round-off  problem  might  be 
encountered.  However,  this  is  to  be  expected  with  either  of  the  two  algorithms. 
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The  second  algorithm  is  a  simple  example  of  a  recursive  mode  of  operation.  The  key  element  in 
any  recursive  procedure  is  the  use  of  the  results  of  the  previous  step  to  aid  in  obtaining  the  desired  result 
for  the  current  step,  and  is  one  of  the  main  features  of  the  Kalman  filter. 

Now  consider  the  problem  of  estimating  the  state  S  of  a  linear,  dynamic,  discrete-time  system  from  the 
measurements  collected  in  a  finite  observation  interval,  namely  Zk  =  {Z0,Zi,...,Zk}.  The  model  assumed 
for  the  system  is  described  in  equations  (2A.4)  and  (2A.5).  The  a-priori  knowledge  about  this  system 
consists  of  the  following  information: 

the  initial  state  S0  is  a  random  vector,  having  known  mean  value  p0  and  covariance  matrix  P0  >  0, 
the  deterministic  input  Uk,  if  present,  is  known, 

the  random  forcing  input  Gk-Vk  is  a  white  noise  process,  having  zero  mean  and  known  covariance  matrix 

Qk>0, 

the  measurement  error  process  Lk-Wk  is  zero-mean  white  noise  having  covariance  matrix  Rk  >  0, 

the  initial  state  So  is  assumed  to  be  uncorrelated  with  the  disturbances  Vk,  Wk, 

the  noise  processes  Vk,  Wk  are  mutually  uncorrelated,  i.e.  E{Vk-WkT}  =  0. 

Two  approaches  can  be  followed  to  derive  the  Kalman  filter  equations.  The  first  assumes  that  the 
processes  Vk,  Wk  and  the  initial  state  S0  are  normally  distributed,  and  the  mean  square  error  criterion  is 
then  used  to  obtain  an  optimum  filter  which  estimates  Sk  as  Sk=  E{Sk|Zk}.  From  the  assumptions  of  a 
linear  model  and  Gaussian  statistics,  a  linear  filter  results.  The  second  approach  is  that  of  linear  minimum 
mean  square  error  (LMSE),  with  no  assumption  on  the  process  distribution  functions.  It  makes  use  of  the 
following  fundamental  features  of  LMSE  estimation: 

the  orthogonality  principle  between  the  estimation  error  and  the  observations, 
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the  concept  of  the  innovation  sequence, 


the  chance  of  deriving  a  recursive  updating  formula  for  the  estimate  when  a  new  measurement  becomes 
available. 

From  the  a-priori  knowledge  about  the  process  S,  not  only  the  initial  value  S0  can  be  estimated  S°  =  po 
with  covariance  P0,  but  also  any  successive  value  can  be  predicted  as  Sk  =  pk  from  the  equation: 

M-k+i  =  flVlhc  +  BkTJk  (1A.20) 

with 

Pk+1  =  Ok-Pk-0)kT  +  Qk  (1A.21) 

These  describe  the  time  evolution  of  the  unconditioned,  a-priori  mean  and  covariance  of  S.  From  them  it 
can  be  observed  that,  even  though  the  initial  guess  is  very  accurate,  a  prediction  error  results  (Pk  >  0)  due 
to  the  noisy  input  Vk.  On  the  other  hand,  if  the  initial  guess  is  poor,  the  prediction  error  remains  severely 
affected  by  it  for  a  time  depending  on  the  transition  matrix  Ok  even  if  Qk  =  0,  i.e.  no  random  input  is 
present. 

One  potential  problem  of  sequential  estimation  is  termed  divergence.  The  problem  occurs  when 
the  magnitude  of  the  state  vector  covariance  matrix  becomes  relatively  small.  This  decrease  in  AP  occurs 
naturally  as  more  observations  are  processed,  since  the  knowledge  of  the  state  vector  increases  with  the 
number  of  observations  processed.  When  AP  becomes  relatively  small,  the  Kalman  gain  becomes 
correspondingly  small,  since  AP  is  a  multiplicative  factor  in  its  calculation.  The  result  is  that  the 
estimator  ignores  new  data  and  does  not  make  significant  improvements  to  the  state  vector.  While  this 
would  seem  to  be  a  desirable  result  of  the  estimation  process,  sometimes  AP  becomes  artificially  small, 
resulting  in  the  filter  disregarding  valid  observations.  To  correct  this  divergence  problem,  process  noise 
is  introduced  to  limit  AP  from  below. 
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We  now  describe  the  particular  implementation  of  our  model.  The  measurement  vector  in 
Cartesian  coordinates  is, 

(1A.22) 


and  our  state  vector  is, 

M  (1A.23) 

s=  x 

y 

jy 

where  x,  y  represent  the  velocity  of  x  and  y,  respectively.  The  state  transition  matrix  is, 

1  T  0  0")  (1A.24) 

0  10  0 
0  0  1  T 

0  0  0  lj 

where  T  is  the  scan  time  of  the  rotating  radar  antenna.  For  simplicity,  we  assume  that  the  scan  time  is 
identical  to  the  time  between  observations.  For  a  target  moving  in  other  than  path  radial  to  the  radar,  the 
observations  will  be  made  at  times  greater  than  or  less  than  the  scan  time,  depending  on  the  target 
velocity  and  path,  and  on  the  antenna  rotation  direction.  We  choose  a  scan  time  often  seconds  as  a 
reasonable  value. 

The  observation  matrix, 

1  0  0  0")  (1A.25) 

0  0  1  0y 

simply  extracts  x  and  y  from  the  state  vector  S  so  that  the  state  may  be  compared  to  the  measurements. 
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The  measurement  noise,  W,  and  the  process  or  plant  noise,  V,  are  2x1  normally  distributed 
vectors  with  known  variance.  V  represents  the  random  acceleration  of  the  target,  and  W  represents 
random  measurement  error  due  primarily  to  thermal  noise  in  the  radar  receiver. 


L  is  a  2x2  identity  matrix,  which  ensures  that  the  measurement  noise  covariance  matrix,  R,  is  of 
appropriate  dimensions.  To  model  random  accelerations  in  the  x  and  y  directions  due  to  model 
uncertainties,  a  4  x  2  matrix,  G,  is  multiplied  by  the  2  x  1  vector  representing  random  accelerations,  V, 
and  introduced  into  the  Kalman  filter’s  state  estimation  calculations.  To  do  this,  G  is  expressed  as 


flT’ 


G  = 


V 


T 

0 

0 


0 

it 
2  1 


(1A.26) 


In  this  way,  when  G  is  multiplied  by  V,  the  random  accelerations,  the  result  is  in  the  same  dimensions  as 
the  state  vector.  Of  course  these  equations  are  simply  the  result  of  Taylor  series  expansion  of  the 
position,  x,  with  respect  to  t,  expressed  in  scalar  form  as. 


x(t)-x0  +v-t  +  \-a-t 2 


(1A.27) 


where  Xo  is  the  initial  position,  v  is  the  velocity,  and  a  is  acceleration. 
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2A.  Appendix  2  Derivation  of  Kalman  Filter  Equations 


2/1. 1  First  Derivation: 

For  simplicity,  designate  the  two  independent  estimates  s*k|k-x  and  zk  by  s\  and  s*2  respectively. 
Designate  s*k|k,  the  optimum  combined  estimate,  by  s*c.  We  desire  to  find  an  optimum  linear  estimate  for 
s  c.  We  can  designate  this  linear  estimate  as 

sc=k^-s\  +k2-s\  (2A.1) 

We  want  this  estimate  s*c  to  be  unbiased,  it  being  assumed  that  s\  and  s*2  are  unbiased.  Designate  s  as 
the  true  value  of  s.  Obtaining  the  mean  of  (2A.  1 ),  it  follows  that  for  the  estimate  to  be  unbiased 

s  =  kx  ■  s  +  k2  ■  s  (2A.2) 


which  becomes, 


1  —  +  k2 


(2A.3) 


Thus  for  the  estimate  to  be  unbiased  we  require, 


k2  =  1  -  kx 


(2A.4) 


Substituting  (2A.4)  into  (2A.1)  yields, 


5*  =  kx  ■  s*  +(!-*,)■  sj 


(2A.5) 


Let  the  variances  of  s*c,  s'i,  and  s*2  be  designated  as  ct2c,  a2),  and  a22.  Then  (2A.5)  yields 


g]  =k\  -cr2  +(1— &,)2  -g2 


2  _2 


(2A.6) 
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To  find  the  ki  that  gives  the  minimum  ct2c,  we  differentiate  (2A.6)  with  respect  to  ki  and  set  the  result  to 
zero,  obtaining, 


2  ■  kx  -<7\  —  2  -  (l  —  Arj )  •  cr|  =  0 


Hence, 


1  2  2 
0\  +  0*2 


Substituting  (2A.8)  into  (2A.5)  yields, 


* 


o*  i 


-Sy  + 


2  2 

cr,  +  o-j 


Rewriting  (2A.9)  yields 


VCTi 


*  A 


'2 ; 


Note  that  substituting  (2A.8)  into  (2A.6)  yields 


O’c  = 


1  1 
— T"  H - r 


v°-i 


'2  J 


(2A.7) 


(2A.8) 


(2A.9) 


(2A.10) 


(2A.11) 


24.2  Second  Derivation: 

The  second  derivation  employs  a  weighted  least-squares  error  estimate  approach.  In  Figure  2A.1 
we  have  two  estimates  zk  and  s*k|k_,  and  desire  to  replace  these  with  a  combined  estimate  s*k|k  that  has  a 
minimum  weighted  least-squares  error.  For  an  arbitrarily  chosen  s*k|k,  there  are  two  errors.  One  is  the 
distance  of  s*k|k  from  zk;  the  other  is  its  distance  of  s*k|k  from  s*k|k,i.  For  the  minimum  least-squares 
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estimate  we  might  minimize  the  sum  of  the  squares  of  the  distances  (errors)  between  the  measurements 
and  the  best  fitting  line  (trajectory)  to  the  measurements[Bro98].  We  would  like  to  does  this  in  a  way  that 
takes  into  consideration  the  accuracies  of  the  measurements.  For  convenience  let  s*k|k.i  be  more  accurate 
than  zk.  In  this  case  it  is  more  important  that  (s*k|k-i  -  s*k|k)2  be  small,  specifically  smaller  than  (zk  -  s*k|k)2. 
This  would  be  achieved  if  in  finding  the  least  sum  of  the  squares  of  each  of  the  two  errors  we  weighted 
the  former  error  by  a  larger  constant  than  the  latter  error.  We  are  thus  obtaining  a  minimization  of  an 
appropriately  weighted  sum  of  the  two  errors  wherein  the  former  receives  a  larger  weighting.  A  logical 
weighting  is  to  weight  each  term  by  1  over  the  accuracy  of  their  respective  estimates  as  the  following 
equation  does: 

£  _  (z*  ~  sk\k  )2  |  (fk\k-\  ~  ^k\h  )  (2A.12) 

VAR{zk)  VAR(s*k^) 

Here  the  error  (s*k|k_i  -  s*k|k)2  is  weighted  by  1  over  the  variance  of  s*k|k-i  and  (zk  -  s*k,k)2  by  1  over  the 
variance  of  zk.  Thus  if  VAR(s*k|k-i)  «  VAR(zk),  then  l/VAR(s*k|k-i)  »  l/VAR(zk)  and  forces  the  error 
(s*k|k-i  -  s*k|k)2  to  be  much  smaller  than  the  error  (zk  -  s*k!k)2  when  minimizing  the  weighted  sum  of  the 
errors  in  (2A.  12).  This  then  forces  s’k|k  to  be  close  to  s*k|k.1(  as  it  should  be.  The  more  accurate  s\|k-i,  the 
closer  s*k|k  is  to  s*k|k-i-  In  (2A.12)  the  two  errors  are  automatically  weighted  according  to  their  importance, 
the  errors  being  divided  by  their  respective  variances.  On  finding  the  s’k|k  that  minimizes  E  of  (2A.12),  a 
weighted  least-squares  estimate  in  stead  of  just  a  least-squares  estimate  is  obtained.  The  s*k|k  that 
minimizes  (2A.12)  is  found  by  differentiating  (2A.12)  with  respect  to  s*k|k  and  setting  the  result  equal  to 
zero,  yielding 


dE_  _  2  ‘  (2k  ■Sqt )  |  2  •  (sk|k_i  ^  (2A.13) 

Kk~  VAR (z*)  f^(4_,) 


Solving  (2A.13)  yields 
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f  zk  _ 1 _  v— — / 

VARisl^ )  VAR(zk )  J  ‘  1  +  1 

VAR(s^ )  VAR(zk) 

Now  let  us  extend  this  derivation  to  the  case  of  vector  measurements  and  states.  As  done  for  the 
one-dimensional  case,  we  weight  the  square  of  the  errors  by  1  over  the  variance  of  the  error.  Thus  our 
error  (or  cost)  function  as  in  (2A.12)  now  denoted  as  J  is  as  follows, 

j  (z,-s;,)2 1  (s;M-s;„)2  <2a.is) 

VAR(Zk)  K47?(S;m) 

Of  course  this  equation  is  only  conceptually  correct;  the  mathematically  correct  equivalent  will  be  given 
shortly. 


Now  we  will  us  (2A.1 5)  to  solve  for  the  new  combined  estimate  S*k|k  that  minimizes  the  weighted 
sum  of  the  squares  of  the  errors.  Conceptually,  this  is  done  just  as  it  was  done  for  the  equivalent  one¬ 
dimensional  (2A.  12).  Specifically,  (2A.  1 5)  is  differentiated  with  respect  to  the  combined  estimate  S*k|k 
with  the  resulting  equation  set  equal  to  zero  in  order  to  solve  for  the  combined  estimate  S*k|k  that 
minimizes  the  weighted  sum  of  the  errors  squared.  Note  that  (2A.  1 5)  is  not  mathematically  correct  when 
operating  on  matrices,  as  will  be  shown. 

When  using  matrix  notation,  the  first  term  on  the  right  of  the  equal  sign  in  (2A.15)  must  be 
written  as 


(Z*-  S^)2 

VAR(Zk) 


«(Z* 


<\k)T-  r^cz.-s;,*) 


Where  the  matrix  Rk  is  the  covariance  matrix  of  Zk,  that  is, 


(2A.16) 
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Rk=COV(Zk) 


(2A.17) 


The  inverse  of  the  covariance  matrix  Rk,  which  is  designated  as  Rk'',  takes  the  place  of  dividing  by  the 
variance  of  Zk  when  dealing  with  matrices.  The  corresponding  correct  form  for  the  second  term  on  the 
right  of  (2A.15)  is 


(S*  - S*  )2 


(2A.18) 


Where  P‘kM  is  the  covariance  matrix  ofS*k,k-r  Substituting  (2A.17)  and  (2A.18)  into  (2A.15)  yields 


J  =  (Zfr -H-Sklk)TR-k'(Zk  -H-S;,)  +  (S;m  -s;|k)  (2A.19) 


Now  we  may  solve  for  the  optimal  S*k|k  as  we  did  in  the  scalar  case,  by  differentiating  the  error 
function  J  with  respect  to  S\Ik,  setting  the  equation  equal  to  zero  and  solving  for  S*k|k.  Differentiation  of  a 
matrix  equation  is  achieved  by  obtaining  the  gradient  of  the  matrix  equation  as  given  by 


Gradient  of  J  s 


dJ 


aj  dj_  dj_ 
ds]  ’  ds2  ’  ’’  dsn 


(2A.20) 


Applying  (2A.20)  to  (2A.19)  yields 


dJ 


dS 


=  2 


k,k 


' (Si,,  - Sl>t_,  )T  ■  P +  2(Z*  - H •  )r  •  R-k'  ■  (-H)  =  0 


(2A.21) 


Where  H  is  the  observation  matrix  which  (in  our  case)  merely  extracts  the  measurement  estimates  from 
S.  Equation  (2A.21)  can  be  rewritten  as 

•  (P*V  +  Hr  .  Rk'  •  H)  =  S5_,  •  p;L,  +  z l .  R,-1 .  H  (2A.22) 

Which  on  taking  the  transpose  of  both  sides  and  using  the  identity  (AB)T=BTAT  yields 
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Off  +Hr  R*’  -H)-S;|t  =P,V_,  -s;,,.,  +  Hr  R*1  -Zk 


(2A.23) 


or, 


s;„  =(p;^  +  Hr  .r,-1  -h)-1  .(p^1.,  -s;m  +nr  -R~k]-zk) 


Applying  the  well-known  matrix  inversion  lemma, 


(p»y +nrR,-!  h)^'  =  p,;_,  -  p,y,  ■  Hr .  (r,  +  h  p(;_,  hVh 


This  can  be  rewritten  as 


(V  +H  'R,-1  -H)-'  =p;m  -K,  H  p;,_, 


where  Kk  is  often  called  the  Kalman  gain  and  is  given  by, 


Kyt  =  P*|*-i  ’Hr  -(R*  +H-PA*|jt_,  -Hr)_1 


Substituting  (2A.26)  into  (2A.24)  yields 


Sk\k  -(P*|t-i  -K/t ' H ' P*|Ar-i )  '  ’(Pffi  -S,M  +Hr  Rk  -Zk) 


'  H  ■  ^k\k-\  +  (P*|4-l  “Kj  •  H  •  P )  •  H  •  R  •  Z 


But  as  will  be  shown  shortly. 


ka=(p;m-k*-h-p;m)-h7,-r-1 


Then  (2A.28)  can  be  rewritten  as 


(2A.24) 


Pff,  ^2A-25) 


(2A.26) 


(2A.27) 


(2A.28) 


(2A.29) 
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(2A.27) 


+Kt  •  Zk  -  Kk  •  H  ■  SJt|A_, 

or 

=S*|*-i+K*  •  (Zk  -H-S^.,) 

Now  we  can  prove  (2A.29).  From  (2A.27)  it  follows  that 

P;M-H7’=Ki.(Ri+H-P;M-Hr) 
Finishing  the  proof,  this  equation  can  be  written  as 

P,V,  Hr  R*1  -K*  -  H.p;M  Hr  R*1  =  Kk 


(2A.28) 


(2A.29) 


(2A.30) 
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3 A.  Appendix  3  Local  to  Global  Coordinate  Transformation 


Figure  3A-1  Target  in  two  Radar  System 


Let  the  radars  Rj  and  R2  be  situated  at  points  Pj  and  P2  on  the  earth’s  surface  as  shown  in  Figure 
3A-1.  An  airplane  A  is  detected  by  the  radar  R,  in  its  local  Cartesian  coordinate  system  P,  with  the 
coordinates  [xi.y^z^sACPi),  and  it  is  necessary  to  find  the  coordinates  [x2,y2,z2]=A(P2)  of  the  airplane  A 
in  the  Cartesian  coordinate  system  P2  of  the  radar  R2.  The  points  P,,  P2  are  given  by  the  geographical 
coordinates  [fiJi]  and  [f2,l2],  where  f,,f2  are  the  geographical  latitudes  and  lbl2  are  the  geographical 
longitudes.  The  xry,  plane  is  tangent  to  the  ellipsoid  at  the  point  P,.  The  Xi  axis  lies  in  the  plane  of  the 
latitude  fj,  the  yj  axis  lies  in  the  plane  of  the  meridian  I,.  The  z,  axis  is  positively  oriented  towards  the 
earth’s  center.  The  x,  axis  and  the  y,  axis  are  oriented  in  the  direction  of  increasing  longitude  and 
latitude,  respectively.  The  axes  of  the  coordinate  system  P2  are  similarly  oriented. 

Next  it  is  desired  to  transform  the  geographical  coordinates  [f,l]  of  a  point  P  into  geocentric 
Cartesian  coordinates  (C;Xc,yc,zc)  =  C,  with  origin  C  at  the  earth’s  center  as  seen  in  Figure  3A-2.  The  zc 
axis  goes  through  the  earth’s  poles,  with  the  positive  orientation  to  the  North  Pole.  The  xc  axis  passes 
through  the  null  and  180th  meridian,  the  positive  orientation  to  the  null  meridian.  The  yc  axis  is 
perpendicular  to  the  xc  axis  and  to  the  zc  axis,  with  positive  orientation  to  the  90th  meridian. 
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Figure  3A-2  Transformation  from  Geographical  to  Geocentric  Coordinates 


The  shape  of  the  Earth  is  a  geoid,  however  for  simplicity,  an  ellipsoid  will  be  used  as  an 
approximation,  see  Figure  3A-3.  Denote  the  major  semi-axis  as  A  =  6378.245  [km]  and  the  minor  semi¬ 
axis  B  =  6356.863  [km].  To  derive  the  transformation  equations  it  is  convenient  to  use  the  parametric 
description  of  the  ellipse.  The  usual  parametric  representation  of  the  ellipse  is 

x(cp)  =  Acos(tp),  z(cp)  =  Bsin(cp)  (3 A.  1) 

In  order  to  relate  the  angle  f  to  x  and  z  we  must  first  compute  the  normal  vector  n  at  [x,z] 


fx  ) 

^-/Isin(^)''' 

=>  w  = 

r-B  cos(^)N 

v  Bcos{cp)  j 

UJ 

K-Asm{(p)j 

Therefore, 


tan(/)  =  --  =  ^tan(». 
z  B 


D 

From  tan(^)  =  — tan(/)  we  can  obtain, 
A 


(3A.3) 
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cos(<p)  = 


•Jl  +  tan2(^) 


,sin(^)  =  tan(^)cos($>) 


(3A.4) 


Figure  3A-3  An  Ellipse  Description  by  the  Slope  of  the  Normal  Vector 


Now  the  transformation  from  the  coordinate  system  Pi  in  to  the  coordinate  system  P2  can  be 
performed,  see  Figure  3A-4.  This  transformation  consists  of  four  partial  transformations,  more  precisely, 
one  translation  and  three  rotations.  Since  the  coordinate  system  Pi  changes  after  each  transformation,  we 
must  distinguish  the  coordinate  systems  using  a  superscript: 


Pj,i  =  0,...,4,;P]°  =P},P/  =P2 


(3A.5) 


3A.1  The  Translation: 
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Figure  3A-4  Coordinate  Transformation  Diagram 


The  values  of  the  displacement  [Ax,Ay,Az]  =  P2(PI)-P](P])  =  A p  are  the  coordinates  of  the  point  P2  in 
the  first  coordinate  system,  denoted  by  P2(Pi).  The  Ap  is  known  in  the  C  system: 

Ap(C)  =  Pj  (C)  -  P\  (C).  Now  we  need  to  transform  Ap(C)  into  Ap(P, ) .  This  transformation  can  be 
computed  as  follows. 

First,  we  must  place  the  origin  Pi"  of  the  coordinate  system  Pi"  at  the  intersection  of  the  equator  and  the 
Greenwich  meridian.  This  transformation  is: 

'0  1  (A  (3A.6) 

,M=  0  0  1. 

c  I"1  0 


70 


To  move  P,"  into  P,',  we  perform  a  rotation  around  the  y"  axis  with  the  angle  1,.  This  transformation  can 
be  described  by  the  rotation  matrix  Ri . 


(Ax') 

r  cos(/, )  0  sin(/])N| 

Ay 

=  R}  x 

Ay 

A  = 

0  1  0 

vAzy 

pj 

vAz; 

p. 

,-sin(/,)  0  cos(/, ) j 

(3A.7) 


To  move  P,'  into  P,,  we  use  a  rotation  around  the  x'  axis  with  the  angle  f,.  This  transformation  is 
described  by  the  matrix  R2: 


(Ax) 

1  1 

(Ax') 

'10  0  ^ 

Ay 

X 

(N 

1! 

Ay 

A" 

0  cos(/, )  sin(/, ) 

1 

A 

pi 

v°  -sin(/i)  cos (/,), 

Now  the  coordinate  system  Pi  is  in  a  general  position.  Summarizing  the  previous  steps,  the  translation 
vector  Ap(Pj)  can  be  expressed  by 


Ay 


=  R2  X  X  M  X 


P2  y  ~  P\  y 
\P2z  ~P\z  J 


(3A.9) 


The  new  coordinates  ACP/)  can  be  obtained  from  the  previous  coordinates  A(P,)  by  using  the  relation 


(x) 

(  Ax) 

A(pi)  = 

y 

=  A(P,)- 

Ay 

vzJ 

Pi 

vAzj 

(3A.10) 


3A.2  First  Rotation: 


71 


We  rotate  the  Pi*  coordinate  system  around  the  x  axes  by  the  angle  f,  in  order  to  make  the  y  axis 
parallel  to  the  earth  axis  zz.  Then,  the  new  coordinates  of  the  airplane  A(P]2)  are  obtained  from  the 
previous  coordinates  A(P/)  using 


AP,2)  = 


(x > 

'l  0  O'] 

y 

-  7?3  X 

y 

,R*  =Rl  ^ 

0  cos(/,)  -  sin(/, ) 

Pi 

KZJ 

p\ 

v°  sinC/i )  cos(/, )  j 

(3A.11) 


3A.3  Second  Rotation: 


We  rotate  the  Pi2  coordinate  system  by  the  angle  A1=12-1i  around  the  y  axis  in  order  to  let  the  two 
x  axes  ofP,2  and  P2  coincide.  The  new  coordinates  A(PI3)  of  the  airplane  are  then  computed  from  the 
previous  coordinates  A(Pi2)  by 


f X ) 

f  cos(A /) 

0 

sin(A/)N| 

y 

X 

-3- 

II 

y 

— 

0 

1 

0 

kzj 

- 

lzv 

Pi 

^-sin(A/) 

0 

cos(A/)v 

(3A.12) 


3A.4  Third  Rotation: 

We  rotate  the  system  Pi3  by  the  angle  f2  around  the  x  axis  to  let  it  coincide  with  the  system  P2. 
Thus,  the  final  coordinates  A(P2)  =  A(Pi4)  of  the  airplane  A  are  obtained  from  the  previous  coordinates 
A(P,3)by 


(X) 

f  1 

0 

0  N 

A{P2)  = 

y 

=  R5  x 

y 

,r5  = 

0 

cos(/2) 

sin(/2 ) 

VZJ 

p/ 

KZ) 

Pi 

,0 

-sin(/2) 

cos(/2); 

(3A.13) 


The  above  steps  can  be  combined  by  concatenation.  By  using  matrices  R]  -  R5  we  get  the  overall 
transformation 
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AP2)  = 


(3A.14) 


fx~) 

f 

( p  _  p  > 
^2^  -rU 

\ 

y 

=  R5  x  x  Rl  x 

A(Pj  )-R2xR]xMx 

P  _  P 
^  i 1^ 

KZJ 

A 

\ 

P  -  P 

\r2 z  r\z  ) 

cy 

3 A.  5  MatLab  Routines 

The  following  two  MatLab  functions  are  (1)  the  S-function  used  for  calling  the  algorithm  that 
performs  the  above  calculations,  and  (2)  the  actual  function  that  performs  the  algorithm. 


function  [sys,xO,str,ts]  =  s_change_view(t,x,u,flag) 
switch  flag, 

%%%%%%%%%%%%%%%%%% 

%  Initialization  % 

%%%%%%%%%%%%%%%%%% 
case  0, 

[sys,xO,str,ts]  =  mdllnitializeSizes; 

%%%%%%%%%% 

%  Update  % 

%%%%%%%%%% 

%case  2, 

%  sys  =  mdlUpdate; 

%%%%%%%%%% 

%  Output  % 

%%%%%%%%%% 
case  3, 

sys  =  mdlOutputs(t,x,u); 

%%%%%%%%%%%%% 

%  Terminate  % 

%%%%%%%%%%%%% 
case  {1,2, 4, 9}, 
sys  =  [];  %  do  nothing 
%%%%%%%%%%%%%%%%%%%% 

%  Unexpected  flags  % 

%%%%%%%%%%%%%%%%%%%% 

otherwise 

error(['unhandled  flag  =  ',num2str(flag)]); 
end 

%end  dsfunc 
% 

%================================================== 

%  mdllnitializeSizes 

%  Return  the  sizes,  initial  conditions,  and  sample  times  for  the  S-function. 

% 

function  [sys,xO,str,ts]  =  mdllnitializeSizes 
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sizes  =  simsizes; 
sizes.NumContStates  =  0; 
sizes.NumDiscStates  =  0; 
sizes.NumOutputs  =  3; 
sizes.Numlnputs  =  - 1 ; 
sizes.DirFeedthrough  =  1; 
sizes.NumSampleTimes  =  1; 
sys  =  simsizes(sizes); 
x0  =  []; 
str  =  []; 
ts  =[-l  0]; 

o/0===========================================:===== 

%  mdlUpdate 

%  Handle  discrete  state  updates,  sample  time  hits,  and  major  time  step 
%  requirements. 

o/0================================================= 

% 

%==================================== 

%  mdlOutputs 

%  Return  the  output  vector  for  the  S-function 
%========================—=============== 

function  sys  =  mdlOutputs(t,x,u) 

%function[A2]=transf(lat  1 ,  Ion  1 ,  A 1  ,lat2,lon2) 

%latl-2,lonl-2  =  geographical  radar  positions  in  (radians) 
%Al=[x;y;z],(P  Imposition  of  airplane  seen  by  radar  at  PI 
%A2=position  of  airplane  seen  by  radar  at  P2 
sys=transf(u(  1  ),u(2),u(3 : 5),u(6),u(7)); 


%sys=transf(lat  1  Jon  1  ,(x,y,z)of  target, lat2,lon2) 

%sys  is  [x;y;z]  position  of  target  seen  by  radar2  at  lat2,lon2 
function[A2]=transf(lat_l,lon_l,A_l,Iat_2,lon_2) 
%latl-2,lonl-2  =  geographical  radar  positions  in  (degrees) 
%Al=[x;y;z],(Pl)=position  of  airplane  seen  by  radar  at  PI 
latl  =deg2rad(lat_l ); 

Ion  1  =deg2rad(lon_l ); 

Iat2=deg2rad(lat_2); 

lon2=deg2rad(lon_2); 

A1=A_1; 

M=[0  1  0; 

00  1; 

-1  0  0]; 

Rl=[cos(lonl)  0  sin(lonl); 

0  10; 

-sin(lonl)  0  cos(lonl)]; 

R2=[l  0  0;0  cos(latl)  sin(latl);0  -sin(latl)  cos(latl)]; 

P 1  c=gg2gc(lat  1  Jon  1 ); 

P2c=gg2gc(lat2,lon2); 
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DPc=P2c-Plc; 

T=R2*Rl*M*DPc; 

a=Al-T; 

Dl=lon2-lonl; 

R4=[cos(Dl)  0  sin(Dl); 

0  1  0; 

-sin(Dl)  0  cos(Dl)]; 

R5=[l  0  0; 

0  cos(lat2)  sin(lat2); 

0  -sin(lat2)  cos(lat2)]; 

A2=R5*R4*R2'*a;  %position  of  airplane  seen  by  radar  at  P2 
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4A.  Appendix  4  Modified  Munkres  Algorithm 


%  Modified  Munkres  Optimal  assignment  algorithm,  rows=tracks,  cols=observations 
function  [x_ret,cov_file_ret,z_file_ret]=munkres(x) 

[nrows,ncols]=size(x); 
k=min(nrows,ncols); 
z_file=randint(nrows,ncols,  80:90); 
z_fi  le=char(z_fi  le); 
cov_file=z_file; 

%  Preliminaries:  if  m>n,  for  each  column  subtract  min(column)  from  that  column 
%  if  n>m,  for  each  row  subtract  min(row)  from  that  row 
if  nrows>ncols, 
for  I=l:ncols, 
x(:,I)=x(:,I)-min(x(:,I)); 
m(I)=fmd(x(:,I)==0); 
end; 

[x_ret,cov_file_ret,z_file_ret]=munk_l(x,cov_file,z_file); 

else; 

for  1=1  :nrows, 
x(I,:)=x(I,:)-min(x(I,:)); 
n(I)=find(x(I,:)==0); 
end 

[x_ret,cov_file_ret,z_file_ret]=munk_l(x,cov_file,z_file); 

end 


% - Step  1 - 

function  [x_ret,cov_file_ret,z_file_ret]:=munk_l(x,cov_file,z_file) 

[j,k]=find(x(:,:)==0); 

for  p=l  :length(j),%  for  all  0's,  if  no  0*  in  row  or  column  of  0,  star  the  0 
if  length(find(z_file(j(p),:)=='*'))==0  &  length(find(z_file(:,k(p))=='*’))==0, 
z_file(j(p),k(p))='*'; 
end 
end 

if  ~exist('x_ret') 

[x_ret,cov_file_ret,z_file_ret]=munk_2(x,cov_file,z_file); 

end 

% . . Step  2 - 

function  [x_ret,cov_file_ret,z_file_ret]=munk_2(x,cov_file,z_file) 

if  exist('x_ret') 
else 

while  length(find(any(z_file=='*')==l))<  min(size(x))&  ~exist(’x_ret'),  %go  until  final  solution 
[j,k]=find(z_file(:, :)=='*'); 
forI=l:length(k), 
cov_file(:,k(I))='c'; 
end 

[x_ret,cov_file_ret,z_file_ret]=munk_3(x,cov_file,z_file); 


76 


end 


if  length(find(any(z_file=='*')==l))==  min(size(x))&  ~exist('x_ret'), 

%x_ret=x; 

x_ret=[]; 

%cov_fi  Ie_ret=cov_fi  le; 
cov_file_ret=[]; 
z_flle_ret=z_file; 
end; 


% - Step  3 - 

function  [x_ret,cov_file_ret,z_file_ret]=munk_3(x,cov_fiIe,z_file) 
if  length(find(any(z_file=- *’)==  1  ))==  min(size(x))&  ~exist('x_ret’) 
x_ret=x; 

cov_file_ret=cov_file; 

z_file_ret=z_file; 

else 


G,k]=find(x(:,:)=0  &  cov_ftle(:,:)~='c'  &  z_file(:, &  z_file(:, :)~=,M); 

%  find  uncovered  zeros 
for  zz=l:length(j), 

if  length(j)>0  &  length(fmd(any(z_file=='*')==l))<  min(size(x)); 
z_file(j(zz),k(zz))-'';%  Pick  uncovered  Z  and  prime  it  to  Z' 
if  length(find(z_file(j(zz),:)=- *'))>=1  ;  %if  there  is  a  0*  in  this  row 
%z_file(j(zz),k(zz))-''; 
b=find(z_file(j(zz), :)=='*'); 
cov_file(:,b)-u';%  uncover  col  of  Z' 
cov_file(j(zz),:)='c';%  cover  row  of  Z' 

elseif  length(find(z_file(j(zz),:)=='*'))==0;%  if  no  Z*  in  in  row  of  Z'  do  step  4 

[x_ret,cov_file_ret,z_file_ret]=munk_4(x,cov_file,z_file) 

end 

if  exist('xjret') 
break 
end 
end 

U,k]=find(x(:,:)=0  &  cov_file(:,:)~-c'  &  z_file(:, &  z_file(:, :)~=p'); 
end 

if  length(find(any(z_file=='*')=l  ))<  min(size(x))&  ~exist('x_ret') 
[x_ret,cov_file_ret,z_file_ret]=munk_5(x,cov_file,z_file); 
end 
end 

% - %  Step  4 - 

function  [x_ret,cov_file_ret,z_file_ret]=munk_4(x,cov_file,z_file)%Step  4 
%z_file_seq=z_file;%  temp  file  to  record  sequence  of  z',  z*,  etc 

if  length(find(any(z_file=- *')==1  ))==  min(size(x))&  ~exist('x_ret') 
x_ret=[];%x; 

cov_file_ret=[]  ;%cov_file; 
z_file_ret=z_file; 
else 


77 


[u,v]=find(z_file(:, &  cov_file(:,:)~-c');%  find  the  uncovered  Z' 
uu=find(z_file(:,v)=='*'); 
if  isempty(uu), 
z_file(u,v)='*'; 
else 

z_file(u,v)='*';%  Star  the  Z’  with  a  Z*  in  its  column 
z_file(uu,v)=' %  Un-star  the  Z*  in  that  column 
while  ~isempty(uu),  %  while  any  Z*  in  col  of  Z',  Z'  in  row  of  Z*,  etc 
w=find(z_file(uu, Find  Z'  in  row  of  un-starred  Z* 
if  ~isempty(vv)%  If  there  is  a  Z'  in  row  of  un-starred  Z* 
z_file(uu,vv)- %  Star  the  Z' 
uu=find(z_fi le( : , vv)=- *');%  Find  Z*  in  col  of  Z' 
if  ~isempty(uu)%  if  there  is  a  Z* 
z_file(uu,vv)=' ';%  Blank  it  out 
end 
else 
uu=[]; 
end 
end 
end 

z_file(find(z_file(:, Star  all  Z'  from  sequence 
z_file(fmd(z_file(:, :)=-''))=' ';%  Erase  all  Primes  from  Z' 
cov_file(find(cov_file(:,:)))— u';  %  Uncover  all  covered  rows/cols 
if  length(find(any(z_file=='*')==l  ))==  min(size(x))&  ~exist('x_ret') 
x_ret=[];%x; 

cov_fi  le_ret=  []  ;%cov_fi  le; 
z_file_ret=z_file; 
else 

[x_ret,cov_file_ret,z_file_ret]=munk_2(x,cov_file,z_file); 

end 

end 

%  Modified  Munkres  Optimal  assignment  algorithm,  rows=tracks,  cols=observations 

% - - - . —  Munkres  Algorithm  Step  5 - 

function  [x_ret,cov_file_ret,z_file_ret]=munk_5(x,cov_file,z_file) 
if  length(find(any(z_file=- *')==1  ))==  min(size(x)) 
x_ret=x; 

cov_file_ret=cov_file; 

z_file_ret=z_file; 

else 

h=min(x(find(cov_file(:,:)~-c'))); 
i=all(cov_file=='c',2);%find  covered  row 
i=find(i==l);%  #  of  covered  rows 

if  length(i)>0  &  length(find(any(z_file=='*')==l))<min(size(x)); 
for  I=l:length(i), 

x(i(I), : )=x(i(I),  :)+h;%  add  h  to  covered  rows 
end 
end 

z=any(cov_file~='c');%find  uncovered  columns 
z=find(z==l); 

if  length(z)>0  &  length(find(any(z_file=='*')==l))<  min(size(x)); 
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for  I=l:Iength(z), 
x(:,z(I))=x(:,z(I))-h; 
end 

%  go  back  to  step  3 

[j,k]=find(x(:,:)==0  &  cov_file(:,:)~='c');  %  find  uncovered  zeros 
end 

[x_ret,cov_file_ret,z_file_ret]=munk_3(x,cov_file,z_file); 

end 
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6.  List  of  Notation  /  Symbols 


B,G,L  matrices  in  dynamic  system  model  that  account  for  measurement  and  model  uncertainties 

H  observation  matrix,  (extracts  measurement  variables  from  state  vector  S) 

d>  state  transition  matrix 

E{.}  expectation  operator 

E{x|z}  conditional  expectation  of  x  given  z 

Kic  Kalman  gain  matrix  at  time  k 

P  covariance  matrix 


P  estimation  covariance  matrix 

a  /\ 

covariance  matrix  of  filtered  and  predicted  estimates,  respectively 

Q>R  process  covariance,  and  measurement  covariance  matrices,  respectively 

Sr  system  state  at  the  k"'  time  instant 

/V 

filtered  estimate  of  Sk 

predicted  estimate  of  Sk 
Z  vector  of  measurements 

v  k  innovation  process  at  time  k 

©*  covariance  matrix  of  the  innovation 
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